/*
$Id: triang3.cc,v 1.5 1996/10/11 16:13:22 roitzsch Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "triang3.h"
#include "utils.h"
#include "numerics.h"
#include "mzibutil.h"
#include "triangtempl.h"
#include "cmdpars.h"
extern CmdPars Cmd;
Vector<PT*> EDG3::p(2);
Vector<PT*> TR3::p(3);
Vector<EDG*> TR3::e(3);
Vector<PATCH*> TR3::f(3);
Vector<PT*> TET3::p(4);
Vector<EDG*> TET3::e(6);
Vector<PATCH*> TET3::f(4);
//-------------------------------------------------------------------------
static const Real oneSixth = 1./6.;
//-------------------------------------------------------------------------
StaticAllocator<TR3> TR3::alloc(ElementsInBlock);
//-------------------------------------------------------------------------
MESH3:: MESH3(const char* inFileName, Bool readMesh)
: MESH(inFileName), ptList(0,25), trList(0,25), edgList(0,25), tetList(0,25),
ptAlloc(ElementsInBlock), edgAlloc(6*ElementsInBlock),
tetAlloc(5*ElementsInBlock)
{
noOfPoints = 0;
noOfEdges = 0;
noOfTriangles = 0;
noOfTetrahedra = 0;
// if (Cmd.isTrue("Circle")) circle = True;
// else circle = False;
// if (Cmd.isTrue("Cylinder")) cylinder = True;
// else cylinder = False;
if (Cmd.isSet("refTetrahedron","BeyStabilDiag"))
{
rightHand = False;
refinementType = BeyStabilDiag;
}
else
{
rightHand = True;
refinementType = ShortDiag;
}
if (readMesh) readTriangulation(inFileName);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
EFMESH3:: EFMESH3(const char* inFileName, Bool /*readMesh*/)
: MESH3(inFileName, False), EFTetAlloc(5*ElementsInBlock)
{
readTriangulation(fileName);
}
//-------------------------------------------------------------------------
MESH3:: ~MESH3()
{
int i;
FORALL( trList,i) delete trList[i];
FORALL( ptList,i) delete ptList[i];
FORALL(edgList,i) delete edgList[i];
FORALL(tetList,i) delete tetList[i];
}
//-------------------------------------------------------------------------
void MESH3:: readTriangulation(const char* fileName0)
{
if (Cmd.isTrue("readHyperFormat")) readHyperFormat(fileName0);
else
if (Cmd.isTrue("readGrapeFormat")) readGrapeFormat(fileName0);
else
if (Cmd.isTrue("readZIBFormat")) readZIBFormat(fileName0);
else
if (Cmd.isTrue("readOldFormat")) readOldZIBFormat(fileName0);
else readZIBFormat(fileName0);
}
//-------------------------------------------------------------------------
void MESH3:: readOldZIBFormat(const char* /*fileName0*/)
{
int ptdim, tetdim, bounddim;
int rc, i, k;
int classA;
int dirTrDim, trNo;
char ch, line[256] ;
PT3 *p;
TR3 *t=0;
TET3 *td;
int maxIndex;
DList<PT3>* ptL = new DList<PT3>; ptList.push(ptL);
DList<EDG3>* edgL = new DList<EDG3>; edgList.push(edgL);
DList<TR3>* trL = new DList<TR3>; trList.push(trL);
DList<TET3>* tetL = new DList<TET3>; tetList.push(tetL);
// BufferedParser parser(fileName, "end", CommentFlag);
Parser parser(fileName, "end", CommentFlag);
if (parser.getLine(line) != parser.ok) missingItem("line",fileName);
if (parser.getLine(line) != parser.ok) missingItem("line",fileName);
rc = sscanf(line,"Dimension:(%d,%d,%d,%d)",
&ptdim,&tetdim,&dirTrDim,&bounddim);
maxIndex = ptdim - 1;
checkMaxIndex(maxIndex);
noOfTetrahedra = tetdim;
Vector<PT3*> ptAdr(0,maxIndex); // for point addresses
FORALL(ptAdr,i) ptAdr[i] = 0;
int ptNo, count=0;
double x, y, z;
for (i=1; i<= ptdim; i++)
{
if (parser.getLine(line) != parser.ok) missingItem("line",fileName);
ch = ' ';
rc = sscanf(line, "%d:(%le,%le,%le),%c%d", &ptNo, &x, &y, &z, &ch,
&classA);
if (rc==0) break;
if (ptNo > ptAdr.high()) pointIndexError(ptNo);
p = ptAlloc.Get(); // new PT3;
ptAdr[ptNo] = p;
p->x = x;
p->y = y;
p->z = z;
p->classA = 0;
if (rc==6) p->classA = classA;
if ((ch=='D')||(ch=='d')) p->boundP = DIRICHLET;
else if ((ch=='I')||(ch=='i')) p->boundP = INTERIOR;
else if ((ch=='N')||(ch=='n')) p->boundP = NEUMANN;
else if ((ch=='C')||(ch=='c')) p->boundP = CAUCHY;
else p->boundP = INTERIOR;
p->node = ++count; // used in Find...
ptList[0]->add(p);
}
noOfPoints=count;
if (parser.findWord("end") != parser.ok) missingItem("end",fileName);
parser.skipRestOfLine();
int pNo1, pNo2 ,pNo3, pNo4, tetNo;
for (i=1; i<= tetdim; i++)
{
classA = 0;
if (parser.getLine(line) != parser.ok) missingItem("line",fileName);
rc = sscanf(line,"%d:(%d,%d,%d,%d),%d", &tetNo,&pNo1,&pNo2,&pNo3,&pNo4,
&classA);
if (rc==0) break;
if ((pNo1>maxIndex)||(pNo2>maxIndex)||(pNo3>maxIndex)||(pNo4>maxIndex))
{
cout << "\n *** ReadTri: bound check error reading tetrahedron\n";
return;
}
td = getTET3(); // new TET3;
td->p1 = ptAdr[pNo1];
td->p2 = ptAdr[pNo2];
td->p3 = ptAdr[pNo3];
td->p4 = ptAdr[pNo4];
td->classA = classA;
SetTPts(td,0);
tetList[0]->add(td);
}
// -- compute edges and faces:
Vector<Stack<EDG3*>*> edsOfPts(noOfPoints);
FORALL(edsOfPts,i) edsOfPts[i] = new Stack<EDG3*>(15);
Vector<Stack<TR3*>*> trsOfPts(noOfPoints);
FORALL(trsOfPts,i) trsOfPts[i] = new Stack<TR3*>(15);
DListIter<TET3> iterTet(tetras());
while ((td=iterTet.all())) completeTetrahedron(td, edsOfPts, trsOfPts);
if (parser.findWord("end") != parser.ok) missingItem("end",fileName);
parser.skipRestOfLine();
PT3 *pt1, *pt2, *pt3;
for (i=1; i<= dirTrDim; i++)
{
if (parser.getLine(line) != parser.ok) missingItem("line",fileName);
ch = ' ';
classA = 0;
rc = sscanf(line,"%d:(%d,%d,%d),%c%d,(%le,%le,%le)",
&trNo,&pNo1,&pNo2,&pNo3,&ch,&classA,&x,&y,&z);
if (rc==0) break;
pt1 = ptAdr[pNo1];
pt2 = ptAdr[pNo2];
pt3 = ptAdr[pNo3];
Stack<TR3*>& trs = *trsOfPts[pt1->node];
FORALL(trs,k)
{
t = trs[k];
if (PinTr(pt1,t) && PinTr(pt2,t) && PinTr(pt3,t)) break;
}
if (!t->onBoundary()) boundaryWarning(pNo1,pNo2,pNo3);
t->classA = classA;
if ((ch == 'd') || (ch == 'D'))
{
t->boundP = DIRICHLET;
(t->e1)->boundP = DIRICHLET;
(t->e2)->boundP = DIRICHLET;
(t->e3)->boundP = DIRICHLET;
(t->e1)->classA = classA;
(t->e2)->classA = classA;
(t->e3)->classA = classA;
}
else if ( (ch == 'c') || (ch == 'C'))
{
t->boundP = CAUCHY;
if ((t->e1)->boundP != DIRICHLET)
{
(t->e1)->boundP = CAUCHY;
(t->e1)->classA = classA;
}
if ((t->e2)->boundP != DIRICHLET)
{
(t->e2)->boundP = CAUCHY;
(t->e2)->classA = classA;
}
if ((t->e3)->boundP != DIRICHLET)
{
(t->e3)->boundP = CAUCHY;
(t->e3)->classA = classA;
}
}
else if ( (ch == 'n') || (ch == 'N'))
{
t->boundP = NEUMANN;
if ((t->e1)->boundP != DIRICHLET)
{
(t->e1)->boundP = NEUMANN;
(t->e1)->classA = classA;
}
if ((t->e2)->boundP != DIRICHLET)
{
(t->e2)->boundP = NEUMANN;
(t->e2)->classA = classA;
}
if ((t->e3)->boundP != DIRICHLET)
{
(t->e3)->boundP = NEUMANN;
(t->e3)->classA = classA;
}
}
}
DListIter<TR3> iterTr(triangles());
while ((t=iterTr.all())) setNeighbsOfTetra(t);
iterTr.reset();
while ((t=iterTr.all()))
if ( !(t->onBoundary()) )
{
trList[0]->remove(t);
delete t;
}
iterTet.reset();
while ((td=iterTet.all())) td->setBoundary();
FORALL(edsOfPts,i) delete edsOfPts[i];
FORALL(trsOfPts,i) delete trsOfPts[i];
if (infoRefinement) Info();
if (Cmd.isTrue("printMesh")) cout << *this;
return;
}
//-------------------------------------------------------------------------
void MESH3:: readGrapeFormat(const char* /*fileName0*/)
{
int i, j, maxIndex, noOfTds;
int count=0;
int classA, pNo1, pNo2 ,pNo3, pNo4;
PT3 *p, *pt1=0, *pt2=0, *pt3=0;
TR3* t=0;
TET3* td;
DList<PT3>* ptL = new DList<PT3>; ptList.push(ptL);
DList<TET3>* tetL = new DList<TET3>; tetList.push(tetL);
DList<TR3>* trL = new DList<TR3>; trList.push(trL);
DList<EDG3>* edgL = new DList<EDG3>; edgList.push(edgL);
BufferedParser parser(fileName, "end", CommentFlag);
parser.getValue(&noOfTds);
parser.getValue(&maxIndex);
checkMaxIndex(maxIndex);
Vector<PT3*> ptAdr(0,maxIndex); // for point addresses
FORALL(ptAdr,i) ptAdr[i] = 0;
for (i=1; i <= maxIndex; i++)
{
p = ptAlloc.Get(); // new PT3;
ptAdr[i] = p;
parser.getValue(&p->x);
parser.getValue(&p->y);
parser.getValue(&p->z);
p->node = ++count; // used in Find...
ptList[0]->add(p);
++noOfPoints;
}
rightHand = False;
for (i=1; i <= noOfTds; i++)
{
parser.getValue(&pNo1);
parser.getValue(&pNo2);
parser.getValue(&pNo3);
parser.getValue(&pNo4);
parser.getValue(&classA);
td = getTET3(); // new TET3;
td->p1 = ptAdr[pNo1];
td->p2 = ptAdr[pNo2];
td->p3 = ptAdr[pNo3];
td->p4 = ptAdr[pNo4];
td->classA = 1;
SetTPts(td,0);
tetList[0]->add(td);
++noOfTetrahedra;
}
// -- compute edges and faces:
Vector<Stack<EDG3*>*> edsOfPts(noOfPoints);
FORALL(edsOfPts,i) edsOfPts[i] = new Stack<EDG3*>(15);
Vector<Stack<TR3*>*> trsOfPts(noOfPoints);
FORALL(trsOfPts,i) trsOfPts[i] = new Stack<TR3*>(15);
DListIter<TET3> iterTet(tetras());
while ((td=iterTet.all())) completeTetrahedron(td, edsOfPts, trsOfPts);
int no, pNo[5];
iterTet.reset();
while ((td=iterTet.all()))
{
parser.getValue(&no); pNo[1]=no;
parser.getValue(&no); pNo[2]=no;
parser.getValue(&no); pNo[3]=no;
parser.getValue(&no); pNo[4]=no;
for (j=1; j<=4; j++)
{
if (pNo[j] != -1) continue;
if ( j == 1)
{ pt1 = td->p2; pt2 = td->p3; pt3 = td->p4; }
else if ( j == 2)
{ pt1 = td->p1; pt2 = td->p3; pt3 = td->p4; }
else if ( j == 3)
{ pt1 = td->p1; pt2 = td->p2; pt3 = td->p4; }
else if ( j == 4)
{ pt1 = td->p1; pt2 = td->p2; pt3 = td->p3; }
Stack<TR3*>& trs = *trsOfPts[pt1->node];
FORALL(trs,i)
{
t = trs[i];
if (PinTr(pt1,t) && PinTr(pt2,t) && PinTr(pt3,t)) break;
}
if (!t->onBoundary()) boundaryWarning(pNo1,pNo2,pNo3);
t->boundP = DIRICHLET; t->classA = 1;
t->e1->boundP = t->e2->boundP = t->e3->boundP = DIRICHLET;
t->e1->classA = t->e2->classA = t->e3->classA = 1;
t->p1->boundP = t->p2->boundP = t->p3->boundP = DIRICHLET;
t->p1->classA = t->p2->classA = t->p3->classA = 1;
}
}
DListIter<TR3> iterTr(triangles());
while ((t=iterTr.all())) setNeighbsOfTetra(t);
iterTr.reset();
while ((t=iterTr.all()))
if ( !(t->onBoundary()) )
{
trList[0]->remove(t);
delete t;
}
iterTet.reset();
while ((td=iterTet.all())) td->setBoundary();
FORALL(edsOfPts,i) delete edsOfPts[i];
FORALL(trsOfPts,i) delete trsOfPts[i];
if (infoRefinement) Info();
if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------
void MESH3:: readZIBFormat(const char* /*fileName0*/)
{
int i, maxIndex;
TR3* t;
TET3* td;
DList<PT3>* ptL = new DList<PT3>; ptList.push(ptL);
DList<TET3>* tetL = new DList<TET3>; tetList.push(tetL);
DList<TR3>* trL = new DList<TR3>; trList.push(trL);
DList<EDG3>* edgL = new DList<EDG3>; edgList.push(edgL);
BufferedParser parser(fileName, "end", CommentFlag);
if (parser.findWord("Poin") != parser.ok) missingItem("points",fileName);
if (parser.findWord("max") != parser.ok) missingItem("maxIndex",fileName);
parser.getValue(&maxIndex);
checkMaxIndex(maxIndex);
Vector<PT3*> ptAdr(0,maxIndex); // for point addresses
FORALL(ptAdr,i) ptAdr[i] = 0;
readPoints(parser, ptAdr);
if (parser.findWord("Elem") != parser.ok) missingItem("tetrahedra",fileName);
readElements(parser, ptAdr);
// -- compute edges and faces:
Vector<Stack<EDG3*>*> edsOfPts(noOfPoints);
FORALL(edsOfPts,i) edsOfPts[i] = new Stack<EDG3*>(15);
Vector<Stack<TR3*>*> trsOfPts(noOfPoints);
FORALL(trsOfPts,i) trsOfPts[i] = new Stack<TR3*>(15);
DListIter<TET3> iterTet(tetras());
while ((td=iterTet.all())) completeTetrahedron(td, edsOfPts, trsOfPts);
if (parser.findWord("Boun") != parser.ok) missingItem("boundary",fileName);
readBoundary(parser, ptAdr, trsOfPts);
DListIter<TR3> iterTr(triangles());
while ((t=iterTr.all())) setNeighbsOfTetra(t);
iterTr.reset();
while ((t=iterTr.all()))
if ( !(t->onBoundary()) )
{
trList[0]->remove(t);
delete t;
}
iterTet.reset();
while ((td=iterTet.all())) td->setBoundary();
FORALL(edsOfPts,i) delete edsOfPts[i];
FORALL(trsOfPts,i) delete trsOfPts[i];
if (infoRefinement) Info();
if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------
void MESH3:: readPoints(BufferedParser& parser, Vector<PT3*>& ptAdr)
{
int ptNo, count=0;
PT3* p;
while (parser.getValue(&ptNo) == parser.ok)
{
if (ptNo > ptAdr.high()) pointIndexError(ptNo);
p = ptAlloc.Get(); // new PT3;
ptAdr[ptNo] = p;
parser.getValue(&p->x);
parser.getValue(&p->y);
parser.getValue(&p->z);
p->node = ++count; // used in Find...
ptList[0]->add(p);
++noOfPoints;
}
}
//-------------------------------------------------------------------------
void MESH3:: readElements(BufferedParser& parser, Vector<PT3*>& ptAdr)
{
int pNo1, pNo2 ,pNo3, pNo4, classA;
TET3 *td;
while (parser.getValue(&pNo1) == parser.ok)
{
parser.getValue(&pNo2);
parser.getValue(&pNo3);
parser.getValue(&pNo4);
td = getTET3(); // new TET3;
td->p1 = ptAdr[pNo1];
td->p2 = ptAdr[pNo2];
td->p3 = ptAdr[pNo3];
td->p4 = ptAdr[pNo4];
parser.getValue(&classA); // int -> char
td->classA = classA;
SetTPts(td,0);
tetList[0]->add(td);
++noOfTetrahedra;
}
}
//-------------------------------------------------------------------------
void MESH3:: readBoundary(BufferedParser& parser, Vector<PT3*>& ptAdr,
Vector<Stack<TR3*>*>& trsOfPts)
{
int i, pNo1, pNo2 ,pNo3, pNo;
PT3 *pt1, *pt2, *pt3, *p;
TR3 *t = 0;
while (parser.getValue(&pNo1) == parser.ok) // triangles on boundaries
{
parser.getValue(&pNo2);
parser.getValue(&pNo3);
pt1 = ptAdr[pNo1];
pt2 = ptAdr[pNo2];
pt3 = ptAdr[pNo3];
Stack<TR3*>& trs = *trsOfPts[pt1->node];
FORALL(trs,i)
{
t = trs[i];
if (PinTr(pt1,t) && PinTr(pt2,t) && PinTr(pt3,t)) break;
}
if (!t->onBoundary()) boundaryWarning(pNo1,pNo2,pNo3);
t->readBC(parser);
}
if (parser.findWord("Poin") == parser.ok) // additional BCs on points
{
while (parser.getValue(&pNo) == parser.ok)
{
p = ptAdr[pNo];
p->readBC(parser);
}
}
}
//-------------------------------------------------------------------------
void MESH3:: readHyperFormat(const char* /*fileName0*/)
{
int i, nPoints, nEdges, nTriangles, nTetrahedrons;
int pointsSection, trianglesSection, tetrahedronsSection;
int pointDataSection, triangleDataSection, tetrahedronDataSection;
TR3* t;
TET3* td;
DList<PT3>* ptL = new DList<PT3>; ptList.push(ptL);
DList<TET3>* tetL = new DList<TET3>; tetList.push(tetL);
DList<TR3>* trL = new DList<TR3>; trList.push(trL);
DList<EDG3>* edgL = new DList<EDG3>; edgList.push(edgL);
BufferedParser parser(fileName, "end", CommentFlag);
if (parser.findWord("nPoints") != parser.ok) missingItem("nPoints",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
parser.getValue(&nPoints);
checkMaxIndex(nPoints);
noOfPoints = nPoints;
if (parser.findWord("nEdges") != parser.ok) missingItem("nEdges",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
parser.getValue(&nEdges);
checkMaxIndex(nEdges);
if (parser.findWord("nTriangles") != parser.ok)
missingItem("nTriangles",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
parser.getValue(&nTriangles);
checkMaxIndex(nTriangles);
if (parser.findWord("nTetrahedrons") != parser.ok)
missingItem("nTetrahedrons",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
parser.getValue(&nTetrahedrons);
//cout << "\n nTetrahedrons = " << nTetrahedrons << "\n";
checkMaxIndex(nTetrahedrons);
noOfTetrahedra = nTetrahedrons;
if (parser.findWord("Points") != parser.ok) missingItem("Points",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(&pointsSection);
//cout << "\n pointsSection = " << pointsSection << "\n";
checkMaxIndex(pointsSection);
if (parser.findWord("Triangles") != parser.ok)
missingItem("Triangles",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(&trianglesSection);
//cout << "\n trianglesSection = " << trianglesSection << "\n";
checkMaxIndex(trianglesSection);
if (parser.findWord("Tetrahedrons") != parser.ok)
missingItem("Tetrahedrons",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(&tetrahedronsSection);
//cout << "\n tetrahedronsSection = " << tetrahedronsSection << "\n";
checkMaxIndex(tetrahedronsSection);
if (parser.findWord("PointData") != parser.ok)
missingItem("PointData",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(&pointDataSection);
//cout << "\n pointDataSection = " << pointDataSection << "\n";
checkMaxIndex(pointDataSection);
if (parser.findWord("TriangleData") != parser.ok) missingItem("TriangleData",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(&triangleDataSection);
//cout << "\n triangleDataSection = " << triangleDataSection << "\n";
checkMaxIndex(triangleDataSection);
if (parser.findWord("TetrahedronData") != parser.ok) missingItem("TetrahedronData",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(&tetrahedronDataSection);
//cout << "\n tetrahedronDataSection = " << tetrahedronDataSection << "\n";
checkMaxIndex(tetrahedronDataSection);
Vector<PT3*> ptAdr(0,nPoints); // for point addresses
FORALL(ptAdr,i) ptAdr[i] = 0;
readHyperPoints(parser, pointsSection, ptAdr);
readHyperElements(parser, tetrahedronsSection, ptAdr);
readHyperClass(parser, tetrahedronDataSection);
// -- compute edges and faces:
Vector<Stack<EDG3*>*> edsOfPts(noOfPoints);
FORALL(edsOfPts,i) edsOfPts[i] = new Stack<EDG3*>(15);
Vector<Stack<TR3*>*> trsOfPts(noOfPoints);
FORALL(trsOfPts,i) trsOfPts[i] = new Stack<TR3*>(15);
DListIter<TET3> iterTet(tetras());
while ((td=iterTet.all())) completeTetrahedron(td, edsOfPts, trsOfPts);
//cout << "\n readHyperBoundary\n";
readHyperBoundary(parser, trianglesSection, triangleDataSection, ptAdr,
trsOfPts);
DListIter<TR3> iterTr(triangles());
while ((t=iterTr.all())) setNeighbsOfTetra(t);
iterTr.reset();
while ((t=iterTr.all()))
if ( !(t->onBoundary()) )
{
trList[0]->remove(t);
delete t;
}
iterTet.reset();
while ((td=iterTet.all())) td->setBoundary();
FORALL(edsOfPts,i) delete edsOfPts[i];
FORALL(trsOfPts,i) delete trsOfPts[i];
if (infoRefinement) Info();
if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------
void MESH3:: readHyperPoints(BufferedParser& parser, int pointsSection, Vector<PT3*>& ptAdr)
{
int count=0;
PT3* p;
int section = 0;
while (section != pointsSection)
{
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(§ion);
}
for (int i=1; i<=noOfPoints; i++)
{
p = ptAlloc.Get(); // new PT3;
ptAdr[i] = p;
parser.getValue(&p->x);
parser.getValue(&p->y);
parser.getValue(&p->z);
p->node = ++count; // used in Find...
// cout << "\n p = " << p->x << " , " << p->y << " , " << p->z << "\n";
ptList[0]->add(p);
}
}
//-------------------------------------------------------------------------
void MESH3:: readHyperElements(BufferedParser& parser, int tetrahedronsSection, Vector<PT3*>& ptAdr)
{
int pNo1, pNo2 ,pNo3, pNo4;
TET3 *td;
int section = 0;
while (section != tetrahedronsSection)
{
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(§ion);
}
for (int i=1; i<=noOfTetrahedra; i++)
{
parser.getValue(&pNo1);
parser.getValue(&pNo2);
parser.getValue(&pNo3);
parser.getValue(&pNo4);
td = getTET3(); // new TET3;
td->p1 = ptAdr[pNo1];
td->p2 = ptAdr[pNo2];
td->p3 = ptAdr[pNo3];
td->p4 = ptAdr[pNo4];
SetTPts(td,0);
// cout << "\n tet = " << pNo1 << ", " << pNo2 << ", " << pNo3 << ", " << pNo4 << "\n";
tetList[0]->add(td);
}
}
//-------------------------------------------------------------------------
void MESH3:: readHyperClass(BufferedParser& parser, int tetrahedronDataSection)
{
int classA;
TET3 *td;
int section = 0;
while (section != tetrahedronDataSection)
{
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(§ion);
}
DListIter<TET3> iterTet(tetras());
while ((td=iterTet.all()))
{
parser.getValue(&classA); // int -> char
td->classA = classA;
//cout << "\n tet->class = " << classA << "\n";
}
}
//-------------------------------------------------------------------------
void MESH3:: readHyperBoundary(BufferedParser& parser, int trianglesSection,
int triangleDataSection,
Vector<PT3*>& ptAdr, Vector<Stack<TR3*>*>& trsOfPts)
{
int i, pNo1, pNo2 ,pNo3, pNo4, pNo5;
PT3 *pt1, *pt2, *pt3;
TR3 *t = 0;
Vector<int> trClass(noOfTriangles), trBound(noOfTriangles);
parser.rewind();
if (parser.findWord("TetrahedronData") != parser.ok) missingItem("TetrahedronData",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
int section = 0;
while (section != triangleDataSection)
{
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(§ion);
}
for (i=1; i<=noOfTriangles; i++) // triangles on boundaries
{
parser.getValue(&pNo1);
parser.getValue(&pNo2);
trBound[i] = pNo1;
trClass[i] = pNo2;
// cout << "\n class, boundP = " << pNo1 << ", " << pNo2 << "\n";
}
parser.rewind();
if (parser.findWord("TetrahedronData") != parser.ok) missingItem("TetrahedronData",fileName);
if (parser.findWord("=") != parser.ok) missingItem("=",fileName);
section = 0;
while (section != trianglesSection)
{
if (parser.findChar('@') != parser.ok) missingItem("@",fileName);
parser.getValue(§ion);
}
for (int j=1; j<=noOfTriangles; j++) // triangles on boundaries
{
parser.getValue(&pNo1);
parser.getValue(&pNo2);
parser.getValue(&pNo3);
parser.getValue(&pNo4);
parser.getValue(&pNo5);
// cout << "\ntriangles = " << pNo1 << pNo2 << pNo3 << pNo4 << pNo5 << "\n";
pt1 = ptAdr[pNo1];
pt2 = ptAdr[pNo2];
pt3 = ptAdr[pNo3];
Stack<TR3*>& trs = *trsOfPts[pt1->node];
FORALL(trs,i)
{
t = trs[i];
if (PinTr(pt1,t) && PinTr(pt2,t) && PinTr(pt3,t)) break;
}
if ( t->onBoundary() )
{
t->classA = trClass[j];
t->boundP = trBound[j];
t->readHyperBC(parser);
}
}
/*
if (parser.findWord("Poin") == parser.ok) // additional BCs on points
{
while (parser.getValue(&pNo) == parser.ok)
{
p = ptAdr[pNo];
p->readBC(parser);
}
}
*/
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool MESH3:: pointOnEdge(PT3* p, EDG3* ed)
{
return (ed->p1 == p || ed->p2 == p);
}
//-------------------------------------------------------------------------
Bool MESH3:: edgeOnTri(EDG3* ed, TR3* t)
{
return (t->e1 == ed || t->e2 == ed || t->e3 == ed);
}
//-------------------------------------------------------------------------
EDG3* MESH3:: FindEdgByPts(PT3 *p1, PT3 *p2, Vector<Stack<EDG3*>*>& edsOfPts)
{
int i;
EDG3 *ed;
Stack<EDG3*>& edges = *edsOfPts[p1->node];
FORALL(edges,i)
{
ed = edges[i];
if (pointOnEdge(p1,ed) && pointOnEdge(p2,ed)) return ed;
}
ed = edgAlloc.Get(); // new EDG3;
ed->p1 = p1;
ed->p2 = p2;
edsOfPts[p1->node]->push(ed);
edsOfPts[p2->node]->push(ed);
edgList[0]->add(ed);
++noOfEdges;
return ed;
}
//-------------------------------------------------------------------------
TR3* MESH3:: FindTrByEds(PT3 *p1, PT3 *p2, PT3 *p3, EDG3 *e1, EDG3 *e2, EDG3 *e3,
Vector<Stack<TR3*>*>& trsOfPts)
{
int i;
TR3 *t;
Stack<TR3*>& trs = *trsOfPts[p1->node];
FORALL(trs,i)
{
t = trs[i];
if (edgeOnTri(e1,t) && edgeOnTri(e2,t) && edgeOnTri(e3,t)) return t;
}
t = new TR3;
noOfTriangles++;
t->e1 = e1;
t->e2 = e2;
t->e3 = e3;
SetTP(t,0);
trsOfPts[p1->node]->push(t);
trsOfPts[p2->node]->push(t);
trsOfPts[p3->node]->push(t);
trList[0]->add(t);
return t;
}
//-------------------------------------------------------------------------
int MESH3:: completeTetrahedron(TET3 *td, Vector<Stack<EDG3*>*>& edsOfPts,
Vector<Stack<TR3*>*>& trsOfPts)
{
PT3 *p1 = td->p1, *p2 = td->p2, *p3 = td->p3, *p4 = td->p4;
EDG3 *e1, *e2, *e3, *e4, *e5, *e6;
TR3 *t1, *t2, *t3, *t4;
td->e1 = FindEdgByPts(p1,p3,edsOfPts);
td->e2 = FindEdgByPts(p2,p3,edsOfPts);
td->e3 = FindEdgByPts(p1,p2,edsOfPts);
td->e4 = FindEdgByPts(p3,p4,edsOfPts);
td->e5 = FindEdgByPts(p2,p4,edsOfPts);
td->e6 = FindEdgByPts(p1,p4,edsOfPts);
e1 = td->e1; e2 = td->e2; e3 = td->e3;
e4 = td->e4; e5 = td->e5; e6 = td->e6;
t1 = FindTrByEds(p2,p3,p4,e2,e4,e5,trsOfPts);
t2 = FindTrByEds(p1,p3,p4,e1,e4,e6,trsOfPts);
t3 = FindTrByEds(p1,p2,p4,e3,e5,e6,trsOfPts);
t4 = FindTrByEds(p1,p2,p3,e1,e2,e3,trsOfPts);
// neighborhood of this tetrahedron
int k = 0;
if ((t1->t31)==0) { t1->t31 = td; k++; }
else if ((t1->t32)==0) { t1->t32 = td; k++; }
if ((t2->t31)==0) { t2->t31 = td; k++; }
else if ((t2->t32)==0) { t2->t32 = td; k++; }
if ((t3->t31)==0) { t3->t31 = td; k++; }
else if ((t3->t32)==0) { t3->t32 = td; k++; }
if ((t4->t31)==0) { t4->t31 = td; k++; }
else if ((t4->t32)==0) { t4->t32 = td; k++; }
if (k!=4)
{
cout << "\n*** Tetra: impossible state (completeTetrahedron), k= "
<< k << "\n";
cout.flush(); abort();
}
return True;
}
//-------------------------------------------------------------------------
void MESH3:: setNeighbsOfTetra(TR3* t)
{
TET3 *tet1 = t->t31, *tet2 = t->t32;
if (tet2 != nil)
{
if (!PinTr(tet2->p1,t))
{
if (tet1) tet2->n1 = tet1;
else { tet2->n1 = t; SetTP(tet2,t,2,3,4); }
}
else
if (!PinTr(tet2->p2,t))
{
if (tet1) tet2->n2 = tet1;
else { tet2->n2 = t; SetTP(tet2,t,1,4,3); }
}
else
if (!PinTr(tet2->p3,t))
{
if (tet1) tet2->n3 = tet1;
else { tet2->n3 = t; SetTP(tet2,t,1,2,4); }
}
else
if (!PinTr(tet2->p4,t))
{
if (tet1) tet2->n4 = tet1;
else { tet2->n4 = t; SetTP(tet2,t,1,3,2); }
}
}
if (tet1 != nil)
{
if (!PinTr(tet1->p1,t))
{
if (tet2) tet1->n1 = tet2;
else { tet1->n1 = t; SetTP(tet1,t,2,3,4); }
}
else
if (!PinTr(tet1->p2,t))
{
if (tet2) tet1->n2 = tet2;
else { tet1->n2 = t; SetTP(tet1,t,1,4,3); }
}
else
if (!PinTr(tet1->p3,t))
{
if (tet2) tet1->n3 = tet2;
else { tet1->n3 = t; SetTP(tet1,t,1,2,4); }
}
else
if (!PinTr(tet1->p4,t))
{
if (tet2) tet1->n4 = tet2;
else { tet1->n4 = t; SetTP(tet1,t,1,3,2); }
}
}
return;
}
//-------------------------------------------------------------------------
void MESH3:: writeTriangulation(const char* /*outFileText*/, const Vector<Real>& u)
const
{
if (Cmd.isTrue("writeGrapeFormat")) writeGrapeFormat (fileName, u);
if (Cmd.isTrue("writeAVSFormat")) writeAVSFormat (fileName, u);
if (Cmd.isTrue("writeHyperFormat")) writeHyperFormat (fileName, u);
if (Cmd.isTrue("writeZIBFormat")) writeZIBFormat (fileName, u);
}
//-------------------------------------------------------------------------
void MESH3:: writeZIBFormat(const char* outFileText, const Vector<Real>& u)const
{
int i;
int node;
TET3 *tet;
TR3 *t;
PT3 *p;
FILE *outFile;
Vector<int> PNodes(noOfPoints);
// saving the node numbers for later reuse
i = 0;
DListIter<PT3> 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<PT3> iterPts(points());
while ((p=iterPts.all())) {
p->node = node+1;
fprintf(outFile, " %d %e %e %e\n", p->node, p->x, p->y, p->z);
node++;
}
fprintf(outFile, "END\n");
fprintf(outFile, "\nElements:\n");
DListIter<TET3> iterTet(tetras());
while ((tet=iterTet.all())) {
fprintf(outFile, "%d %d %d %d %d\n",
(tet->p1)->node, (tet->p2)->node,
(tet->p3)->node, (tet->p4)->node,
tet->classA );
}
fprintf(outFile, "END\n");
fprintf(outFile, "\nBoundary: #faces\n");
DListIter<TR3> iterTrs(triangles());
while ((t=iterTrs.all())) {
if (!t->onBoundary()) continue;
if (t->isDirichlet())
{
if (t->isoP == CIRCULAR)
fprintf(outFile, "%d %d %d D %d S\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
else if (t->isoP == CYLINDRICAL)
fprintf(outFile, "%d %d %d D %d C\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
else fprintf(outFile, "%d %d %d D %d P\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
}
else if (t->isNeumann())
{
if (t->isoP == CIRCULAR)
fprintf(outFile, "%d %d %d N %d S\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
else if (t->isoP == CYLINDRICAL)
fprintf(outFile, "%d %d %d N %d C\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
else fprintf(outFile, "%d %d %d N %d P\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
}
else if (t->isCauchy())
{
if (t->isoP == CIRCULAR)
fprintf(outFile, "%d %d %d C %d S\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
else if (t->isoP == CYLINDRICAL)
fprintf(outFile, "%d %d %d C %d C\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
else fprintf(outFile, "%d %d %d C %d P\n",
(t->p1)->node, (t->p2)->node,(t->p3)->node,t->classA);
}
else
fprintf(outFile, "%d %d %d I \n",
(t->p1)->node, (t->p2)->node,(t->p3)->node );
}
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 ZIB-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 MESH3:: writeHyperFormat(const char* outFileText, const Vector<Real>& u)
const
{
FILE *outFile;
int i;
TET3 *tet;
TR3 *t;
EDG3 *ed;
PT3 *p, *p1, *p2, *p3, *p4;
PATCH *n1, *n2, *n3, *n4;
Vector<int> PNodes(noOfPoints), ENodes(noOfEdges),
TrNodes(noOfTriangles), TetNodes(noOfTetrahedra);
Vector<char> TetMark(noOfTetrahedra);
// saving the node numbers for later reuse and temorary node numbering
i = 0;
DListIter<PT3> PIter(points());
while ((p=PIter.all())) {
PNodes[++i] = p->node;
p->node = i;
}
i = 0;
DListIter<EDG3> EIter(edges());
while ((ed=EIter.all())) {
ENodes[++i] = ed->node;
ed->node = i;
}
i = 0;
DListIter<TR3> TrIter(triangles());
while ((t=TrIter.all())) {
TrNodes[++i] = t->node;
t->node = i;
}
i = 0;
DListIter<TET3> TetIter(tetras());
while ((tet=TetIter.all())) {
TetNodes[++i] = tet->node;
tet->node = i;
TetMark[i] = tet->mark;
tet->mark = 0;
}
char* outText = new char[strlen(outFileText)+10];
strcpy(outText, outFileText);
strcpy(strchr(outText,'.'),".hyper");
outFile = fopen(outText, "w");
if (outFile==0)
{
cout << "\n*** writeHyper: could not open file " << outText
<< "\n";
cout.flush(); abort();
}
fprintf(outFile, "# HyperMesh ASCII 0.2 \n");
fprintf(outFile, "# \n");
fprintf(outFile, "# \n");
fprintf(outFile, "nPoints = %d \n", noOfPoints);
fprintf(outFile, "nEdges = %d \n", noOfEdges);
fprintf(outFile, "nTriangles = %d \n", noOfTriangles);
fprintf(outFile, "nTetrahedrons = %d \n", noOfTetrahedra);
fprintf(outFile, "# \n");
fprintf(outFile, "Points { double[3] } = @1 \n");
fprintf(outFile, "Edges { int[2] } = @2 \n");
fprintf(outFile, "Triangles { int[5] } = @3 \n");
fprintf(outFile, "Tetrahedrons { int[4] } = @4 \n");
fprintf(outFile, "TriangleEncoding = points tetrahedrons \n");
fprintf(outFile, "TetrahedronEncoding = points \n");
fprintf(outFile, "PointData { double } = @5 \n");
fprintf(outFile, "TriangleData { byte[2] } = @6 \n");
fprintf(outFile, "TetrahedronData { byte[2] } = @7 \n");
fprintf(outFile, "# \n");
fprintf(outFile, "#End of Header. Data follows.\n");
fprintf(outFile, "# \n");
fprintf(outFile, "\n@1 \n");
DListIter<PT3> iterPts(points()); // coordinates of the points
while ((p=iterPts.all())) {
fprintf(outFile, "%e %e %e\n", p->x, p->y, p->z);
}
fprintf(outFile, "\n@2 \n");
DListIter<EDG3> iterEds(edges());
while ((ed=iterEds.all())) {
fprintf(outFile, "%d %d\n", (ed->p1)->node, (ed->p2)->node);
}
fprintf(outFile, "\n@3 \n");
int trsMark[1+4];
Vector<char> trBoundP(noOfTriangles), trClassA(noOfTriangles);
DListIter<TET3> iterTets(tetras());
TET3 *nbTet = 0;
int trsM;
i = 0;
while ((tet=iterTets.all())) {
tet->getMarkedTriangles(trsMark);
p1=tet->p1; p2=tet->p2; p3=tet->p3; p4=tet->p4;
n1=tet->n1; n2=tet->n2; n3=tet->n3; n4=tet->n4;
if ( !(tet->orient()) ) {
p1 = tet->p2; p2 = tet->p1;
n1 = tet->n2; n2 = tet->n1;
trsM =trsMark[1]; trsMark[1] = trsMark[2]; trsMark[2] = trsM;
}
if (!trsMark[1])
{
fprintf(outFile, "%d %d %d ", p3->node, p2->node, p4->node);
nbTet = n1->castToTET3();
if ( nbTet )
{
fprintf(outFile, " %d %d", tet->node, nbTet->node );
nbTet->markTriangles(tet);
trBoundP[++i] = INTERIOR; trClassA[i] = 1;
}
else
{
fprintf(outFile, " %d 0", tet->node);
t = n1->castToTR3();
trBoundP[++i] = t->boundP; trClassA[i] = t->classA;
}
fprintf(outFile, "\n");
}
if (!trsMark[2])
{
fprintf(outFile, "%d %d %d ", p1->node, p3->node, p4->node);
nbTet = n2->castToTET3();
if ( nbTet )
{
fprintf(outFile, " %d %d", tet->node, nbTet->node);
nbTet->markTriangles(tet);
trBoundP[++i] = INTERIOR; trClassA[i] = 1;
}
else
{
fprintf(outFile, " %d 0", tet->node);
t = n2->castToTR3();
trBoundP[++i] = t->boundP; trClassA[i] = t->classA;
}
fprintf(outFile, "\n");
}
if (!trsMark[3])
{
fprintf(outFile, "%d %d %d ",
p2->node, p1->node, p4->node);
nbTet = n3->castToTET3();
if ( nbTet )
{
fprintf(outFile, " %d %d", tet->node, nbTet->node);
nbTet->markTriangles(tet);
trBoundP[++i] = INTERIOR; trClassA[i] = 1;
}
else
{
fprintf(outFile, " %d 0", tet->node);
t = n3->castToTR3();
trBoundP[++i] = t->boundP; trClassA[i] = t->classA;
}
fprintf(outFile, "\n");
}
if (!trsMark[4])
{
fprintf(outFile, "%d %d %d ",
p1->node, p2->node, p3->node);
nbTet = n4->castToTET3();
if ( nbTet )
{
fprintf(outFile, " %d %d", tet->node, nbTet->node);
nbTet->markTriangles(tet);
trBoundP[++i] = INTERIOR; trClassA[i] = 1;
}
else
{
fprintf(outFile, " %d 0", tet->node);
t = n4->castToTR3();
trBoundP[++i] = t->boundP; trClassA[i] = t->classA;
}
fprintf(outFile, "\n");
}
}
fprintf(outFile, "\n@4 \n");
iterTets.reset();
while ((tet=iterTets.all()))
{
p1=tet->p1; p2=tet->p2; p3=tet->p3; p4=tet->p4;
if ( !(tet->orient()) ) { p1 = tet->p2; p2 = tet->p1; }
fprintf(outFile, "%d %d %d %d\n", p1->node, p2->node, p3->node, p4->node);
}
fprintf(outFile, "\n@5 \n"); // data at points
iterPts.reset();
while ((p=iterPts.all())) fprintf(outFile, "%e\n", u[PNodes[p->node]]);
fprintf(outFile, "\n@6 \n"); // data at triangles
for (i=1; i<= noOfTriangles; i++)
{
fprintf(outFile, "%d %d \n", trBoundP[i], trClassA[i]);
}
fprintf(outFile, "\n@7 \n"); // data at tetrahedra
iterTets.reset();
while ((tet=iterTets.all()))
fprintf(outFile, "%d\n", tet->classA);
fprintf(outFile, "\n#End of file.\n"); // data at tetrahedra
cout << "\n Hyper-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];
i = 0;
EIter.reset();
while ((ed=EIter.all())) ed->node = ENodes[++i];
i = 0;
TrIter.reset();
while ((t=TrIter.all())) t->node = TrNodes[++i];
i = 0;
TetIter.reset();
while ((tet=TetIter.all()))
{
tet->node = TetNodes[++i];
tet->mark = TetMark[i];
}
return;
}
//-------------------------------------------------------------------------
void MESH3:: writeGrapeFormat(const char* outFileText, const Vector<Real>& u)
const
{
FILE *outFile;
int i;
int node;
int k[4];
TET3 *tet;
PT3 *p;
Vector<int> PNodes(noOfPoints), TetNodes(noOfTetrahedra);
i = 0;
DListIter<PT3> PIter(points());
while ((p=PIter.all())) PNodes[++i] = p->node;
i = 0;
DListIter<TET3> TIter(tetras());
while ((tet=TIter.all())) TetNodes[++i] = tet->node;
char* outText = new char[strlen(outFileText)+10];
strcpy(outText, outFileText);
strcpy(strchr(outText,'.'),".grape");
outFile = fopen(outText, "w");
if (outFile==0)
{
cout << "\n*** grapeTriangulation: could not open file " << outText
<< "\n";
cout.flush(); abort();
}
fprintf(outFile, "%d\n", noOfTetrahedra);
fprintf(outFile, "%d\n", noOfPoints);
node = 0;
DListIter<PT3> iterPts(points());
while ((p=iterPts.all())) {
p->node = node+1;
fprintf(outFile, "%e %e %e\n", p->x, p->y, p->z);
node++;
}
node = 0;
DListIter<TET3> iterTet(tetras());
while ((tet=iterTet.all())) {
tet->node = node+1;
if (tet->orient())
fprintf(outFile, "%d %d %d %d %d\n",
(tet->p1)->node, (tet->p2)->node,
(tet->p4)->node, (tet->p3)->node,
tet->classA );
else
fprintf(outFile, "%d %d %d %d %d\n",
(tet->p1)->node, (tet->p2)->node,
(tet->p3)->node, (tet->p4)->node,
tet->classA );
node++;
}
iterTet.reset();
while ((tet=iterTet.all()))
{
TET3* n;
n = (tet->n1)->castToTET3();
if ( n == nil ) k[0] = -1;
else k[0] = n->node;
n = (tet->n2)->castToTET3();
if ( n == nil ) k[1] = -1;
else k[1] = n->node;
if (tet->orient())
{
n = (tet->n3)->castToTET3();
if ( n == nil ) k[3] = -1;
else k[3] = n->node;
n = (tet->n4)->castToTET3();
if ( n == nil ) k[2] = -1;
else k[2] = n->node;
}
else
{
n = (tet->n3)->castToTET3();
if ( n == nil ) k[2] = -1;
else k[2] = n->node;
n = (tet->n4)->castToTET3();
if ( n == nil ) k[3] = -1;
else k[3] = n->node;
}
fprintf(outFile, "%d %d %d %d\n", k[0], k[1], k[2], k[3]);
}
fprintf(outFile, "1\n");
i = 0;
iterPts.reset();
while ((p=iterPts.all())) fprintf(outFile, "%e\n", u[PNodes[++i]]);
cout << "\n GRAPE-Data ok.; written to " << outText << "\n";
delete outText;
fclose(outFile);
i = 0;
PIter.reset();
while ((p=PIter.all())) p->node = PNodes[++i];
i = 0;
TIter.reset();
while ((tet=TIter.all())) tet->node = TetNodes[++i];
return;
}
//-------------------------------------------------------------------------
void MESH3:: writeAVSFormat(const char* outFileText, const Vector<Real>& u)
const
{
int i;
int node, noOfPValues, noOfTetValues, noOfModValues;
int noOfComponents, noOfElems;
TET3 *tet;
PT3 *p;
FILE *outFile;
Vector<int> PNodes(noOfPoints);
i = 0;
DListIter<PT3> PIter(points());
while ((p=PIter.all())) PNodes[++i] = p->node;
char* outText = new char[strlen(outFileText)+10];
strcpy(outText, outFileText);
strcpy(strchr(outText,'.'),".AVS");
outFile = fopen(outText, "w");
if (outFile==0)
{
cout << "\n*** AVSTriangulation: could not open file " << outText
<< "\n";
cout.flush(); abort();
}
noOfPValues = 1;
noOfTetValues = 1;
noOfModValues = 0;
fprintf(outFile, "%d %d %d %d %d\n",
noOfPoints,noOfTetrahedra,noOfPValues,noOfTetValues, noOfModValues);
node = 0;
DListIter<PT3> iterPts(points());
while ((p=iterPts.all())) {
p->node = node;
fprintf(outFile, "%d %e %e %e\n", p->node, p->x, p->y, p->z);
node++;
}
node = 0;
DListIter<TET3> iterTet(tetras());
while ((tet=iterTet.all())) {
fprintf(outFile, "%d %d tet %d %d %d %d\n",
node, tet->classA,
(tet->p1)->node, (tet->p2)->node,
(tet->p3)->node, (tet->p4)->node
);
node++;
}
noOfComponents = 1;
noOfElems = 1;
fprintf(outFile, "%d %d\n", noOfComponents,noOfElems);
fprintf(outFile, "Temperature, Degree\n");
i = 0;
iterPts.reset();
while ((p=iterPts.all())) {
fprintf(outFile, "%d %e\n", i, u[PNodes[i+1]]);
i++;
}
fprintf(outFile, "%d %d\n", noOfComponents,noOfElems);
fprintf(outFile, "Material, ClassA\n");
i = 0;
iterTet.reset();
while ((tet=iterTet.all())) {
fprintf(outFile, "%d %d\n", i, int(tet->classA));
i++;
}
cout << "\n AVS-Data ok.; written to " << outText << "\n";
delete outText;
fclose(outFile);
i = 0;
PIter.reset();
while ((p=PIter.all())) p->node = PNodes[++i];
}
//-------------------------------------------------------------------------
void MESH3:: Refine()
{
TET3* tet;
int noPt[3], noEdg[3], noTr[3], noTet[3];
Timer timer, accTimer;
noPt[0] = noOfPoints;
noEdg[0] = noOfEdges;
noTr[0] = noOfTriangles;
noTet[0] = noOfTetrahedra;
RemoveGreen(); // delete all green tetras
noPt[1] = noOfPoints;
noTr[1] = noOfTriangles;
noEdg[1] = noOfEdges;
noTet[1] = noOfTetrahedra;
DList<TET3>* tetL = new DList<TET3>;
tetList.push(tetL);
moreRedTd = 0;
DListIter<TET3> iter(tetras());
while ((tet=iter.all())) RefineTetra(tet);
// Algorithm: For all tetrahedra:
// (1) look if there are 3 refined edges,
// two of these on different faces,
// and refine it red this case,
// (2) if no neighbor tetrahedron exists, "father neighbor"
// tetrahedron needs to be red refined (green refined at
// the last level).
// This routine does not use recursion.
if (infoRefinement >= 2) cout << "\n Refine: moreRedTd =";
moreRedTd = 1;
while (moreRedTd)
{
moreRedTd = 0;
iter.reset();
while ((tet=iter.all())) CloseRedTet(tet);
if (infoRefinement >= 2) cout << " " << moreRedTd;
if (moreRedTd == 0) break;
iter.reset();
while ((tet=iter.all())) RefineTetra(tet);
}
noPt[2] = noOfPoints;
noTr[2] = noOfTriangles;
noEdg[2] = noOfEdges;
noTet[2] = noOfTetrahedra;
MakeGreen();
// now: tetrahedra with a boundary face are marked when generated
// iter.reset();
// while (tet=iter.all()) tet->setBoundary();
++refStep;
if (Cmd.isTrue("triCheck"))
consistencyCheck();
if (infoRefinement) refinementInfo (noPt, noEdg, noTr, noTet);
timeInfo(timer, accTimer);
if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------
int MESH3:: CloseRedTet(TET3 *tet)
{
TR3 *t;
TR3 *trs[1+4];
EDG3 *ed, *edSon1, *edSon2;
EDG3 *eds[1+6];
int i, noOfRefEds;
// request to refine a further tetrahedron tet
if ( (tet->mark) == True)
{
moreRedTd++;
return True;
}
noOfRefEds = TestEdsInTd(tet);
if ( noOfRefEds < 3) return True;
if ( noOfRefEds > 3)
{
tet->mark = True;
moreRedTd++;
return True;
}
tet->mark = True;
tet->getTriangles(trs);
for (i=1; i<=4; i++)
if (TestFace(trs[i]) == 3) {tet->mark = False; break;}
if ( tet->mark == True)
{
moreRedTd++;
tet->returnTriangles(trs);
return True;
}
tet->getEdges(eds);
for (i=1; i<=6; i++)
{
ed = eds[i];
edSon1 = ed->firstSon;
if ( (edSon1) != 0)
{
edSon2 = edSon1->next;
if ( (edSon1->firstSon != 0) ||
(edSon2->firstSon != 0) )
{
tet->mark = True;
moreRedTd++;
tet->returnTriangles(trs);
return True;
}
}
}
TET3 *nbTet;
for (i=1; i<=4; i++)
{
t = trs[i]; nbTet = t->t32;
if (nbTet == 0) continue;
if (nbTet->refined())
{
if (testSons(t,nbTet) > 0)
{
tet->mark = True;
moreRedTd++;
tet->returnTriangles(trs);
return True;
}
}
}
tet->returnTriangles(trs);
return True;
}
//-------------------------------------------------------------------------
void MESH3:: RemoveGreen()
{
RelGreen3Td();
RelGreen2ATd();
RelGreen2BTd();
RelGreen1Td();
RelGreen();
RelGreenIrr();
return;
}
//-------------------------------------------------------------------------
int MESH3:: MakeGreen()
{
TET3* tet;
moreGreen1Td = 0;
moreGreen2ATd = 0;
moreGreen2BTd = 0;
moreGreen3Td = 0;
DListIter<TET3> iter(tetras());
while ((tet=iter.all())) MarkGreenTd(tet);
if (infoRefinement >= 2)
cout << "\n GREEN-I -Closure: moreGreen-I-Td = "
<< moreGreen1Td<< "\n";
iter.reset();
while ((tet=iter.all())) {
if (!MkGreen1Td(tet)) {
cout << "*** Tri: failed making GREEN-I tetrahedra (CloseRef)\n";
cout.flush(); abort(); //return False;
}
}
if (infoRefinement >= 2)
cout << " GREEN-II -Closure: moreGreen-IIA-Td = " << moreGreen2ATd;
iter.reset();
while ((tet=iter.all())) {
if (!MkGreen2ATd(tet)) {
cout << "Tri: failed making GREEN-II tetrahedra (CloseRef)\n";
cout.flush(); abort(); // return False;
}
}
if (infoRefinement >= 2)
cout << " moreGreen-IIB-Td = " << moreGreen2BTd << "\n";
iter.reset();
while ((tet=iter.all())) {
if (!MkGreen2BTd(tet)) {
cout << "Tri: failed making GREEN-II tetrahedra (CloseRef)\n";
cout.flush(); abort(); // return False;
}
}
if (infoRefinement >= 2)
cout << " GREEN-III-Closure: moreGreen-III-Td = " << moreGreen3Td <<"\n";
iter.reset();
while ((tet=iter.all())) {
if (!MkGreen3Td(tet)) {
cout << "Tri: failed making GREEN-III tetrahedra (CloseRef)\n";
return False;
}
}
return True;
}
//--------------------------------------------------------------------------------
int MESH3:: RefineTetra(TET3 *td)
{
// internal routine to refine a tetrahedron, loTrues for
if (td->mark) td->mark = False; /* return if not marked */
else return True;
if ((td->firstSon)!=0)
{
cout << "\n*** RefineTetra:??????? marked but already refined \n";
return True; // return if already marked
}
moreRedTd++;
if (refinementType == ShortDiag)
{
if ( !(RefineTdShortDiag(td)) ) {
cout << "\n*** RefineTetra: tetrahedron not refined \n";
cout.flush(); abort();
}
}
else if (refinementType == BeyStabilDiag)
{
if ( !(RefineTdBey(td)) ) {
cout << "\n*** RefineTetra: tetrahedron not refined \n";
cout.flush(); abort();
}
}
else {
cout << "\n*** RefineTetra: unknown refinement type requested\n";
cout.flush(); abort();
}
return True;
}
//-------------------------------------------------------------------------
int MESH3:: RefineTdShortDiag(TET3 *td)
{
// Refinement of a tetrahedron
// Method of Zang (Select the shortest diagonal in the octahedron)
TET3 *td1, *td2, *td3, *td4, *td5, *td6, *td7, *td8;
TR3 *tr, *tr1, *tr2, *tr3, *tr4,
*it1, *it2, *it3, *it4, *it5, *it6, *it7, *it8,
*t11=0, *t12=0, *t13=0, *t14=0,
*t21=0, *t22=0, *t23=0, *t24=0,
*t31=0, *t32=0, *t33=0, *t34=0,
*t41=0, *t42=0, *t43=0, *t44=0,
*t1, *t2, *t3, *t4;
TR3 *trs[1+4];
PT3 *p1 = td->p1,*p2 = td->p2, *p3 = td->p3, *p4 = td->p4,
*p5, *p6, *p7, *p8, *p9, *p10;
EDG3 *e1 = td->e1, *e2 = td->e2, *e3 = td->e3,
*e4 = td->e4, *e5 = td->e5, *e6 = td->e6,
*ie1, *ed;
PATCH *n1=td->n1, *n2=td->n2, *n3=td->n3, *n4=td->n4;
Real len59, len78, len610;
int k, diag59=False, diag78=False, diag610=False;
int depth;
td->getTriangles(trs);
t1=trs[1]; t2=trs[2]; t3=trs[3]; t4=trs[4];
// Get 8 new tetrahedra, 8 new triangles, 1 new edges
// return False if not enough space
td1 = getTET3(); // new TET3;
td2 = getTET3();
td3 = getTET3();
td4 = getTET3();
td5 = getTET3();
td6 = getTET3();
td7 = getTET3();
td8 = getTET3();
noOfTetrahedra += 7;
it1 = new TR3; // new TR3;
it2 = new TR3;
it3 = new TR3;
it4 = new TR3;
it5 = new TR3;
it6 = new TR3;
it7 = new TR3;
it8 = new TR3;
noOfTriangles += 8;
ie1 = edgAlloc.Get(); // new EDG3;
noOfEdges += 1;
if ( (td8==0)||(it8==0)||(ie1==0) )
{
cout << " no memory in RefineTd \n";
return False;
}
depth = (td->depth) + 1;
ie1->depth = depth;
if (ie1->depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[ie1->depth]->add(ie1);
it1->depth = depth;
it2->depth = depth;
it3->depth = depth;
it4->depth = depth;
it5->depth = depth;
it6->depth = depth;
it7->depth = depth;
it8->depth = depth;
// Set the "son" triangles tik, i=1,2,3,4; k=1,2,3,4;
// each triangle is refined
// Triangle 1
if (!RefineTr(t1))
{
cout << " failed in RefineTr called by RefineTd \n";
return False;
}
setTypeOfInnerEdges(t1);
tr1 = t1->firstSon;
if (tr1==0) return False; // No space for Refinement
// Triangle 2
if (!RefineTr(t2))
{
cout << " failed in RefineTr called by RefineTd \n";
return False;
}
tr2 = t2->firstSon;
setTypeOfInnerEdges(t2);
if (tr2==0) return False; // No space for Refinement
// Triangle 3
if (!RefineTr(t3))
{
cout << " failed in RefineTr called by RefineTd \n";
return False;
}
setTypeOfInnerEdges(t3);
tr3 = t3->firstSon;
if (tr3==0) return False; /* No space for Refinement */
// Triangle 4
if (!RefineTr(t4))
{
cout << " failed in RefineTr called by RefineTd \n";
return False;
}
setTypeOfInnerEdges(t4);
tr4 = t4->firstSon;
if (tr4==0) return False; /* No space for Refinement */
p5 = e3->pm;
p6 = e5->pm;
p7 = e2->pm;
p8 = e6->pm;
p9 = e4->pm;
p10= e1->pm;
// one new edge in the octahedron
len59 = Length(p5,p9);
len78 = Length(p7,p8);
len610 = Length(p6,p10);
if ( (len59 <= len78) && (len59 <= len610) )
{
ie1->p1 = p5;
ie1->p2 = p9;
diag59 = True;
}
else if ( (len78 <= len59) && (len78 <= len610) )
{
ie1->p1 = p7;
ie1->p2 = p8;
diag78 = True;
}
else
{
ie1->p1 = p6;
ie1->p2 = p10;
diag610 = True;
}
ie1->pm = 0;
// ie1->bound = 0;
// new faces
tr = tr1;
for (k = 0; k<4; k++)
{
if (PinTr(p2,tr)) t12 = tr;
else if (PinTr(p3,tr)) t13 = tr;
else if (PinTr(p4,tr)) t14 = tr;
else
{
t11 = tr;
ed = EdinTr(tr,p6,p7);
it2->e1 = ed;
if(diag78) it5->e1 = ed;
else if(diag610) it7->e1 = ed;
ed = EdinTr(tr,p7,p9);
it3->e1 = ed;
if(diag59||diag78) it7->e1 = ed;
ed = EdinTr(tr,p6,p9);
it4->e1 = ed;
if(diag59) it5->e1 = ed;
else if(diag610) it8->e2 = ed;
}
tr = tr->next;
}
tr = tr2;
for (k = 0; k<4; k++)
{
if (PinTr(p1,tr)) t21 = tr;
else if (PinTr(p3,tr)) t23 = tr;
else if (PinTr(p4,tr)) t24 = tr;
else
{
t22 = tr;
ed = EdinTr(tr,p8,p10);
it1->e1 = ed;
if(diag78) it8->e1 = ed;
else if(diag610) it6->e1 = ed;
ed = EdinTr(tr,p8,p9);
it4->e2 = ed;
if(diag59) it6->e1 = ed;
else if(diag78) it7->e2 = ed;
ed = EdinTr(tr,p9,p10);
it3->e2 = ed;
if(diag59||diag610) it8->e1 = ed;
}
tr = tr->next;
}
tr = tr3;
for (k = 0; k<4; k++)
{
if (PinTr(p1,tr)) t31 = tr;
else if (PinTr(p2,tr)) t32 = tr;
else if (PinTr(p4,tr)) t34 = tr;
else
{
t33 = tr;
ed = EdinTr(tr,p5,p8);
it1->e2 = ed;
if(diag59||diag78) it6->e2 = ed;
ed = EdinTr(tr,p6,p8);
it4->e3 = ed;
if(diag78) it5->e2 = ed;
else if(diag610) it6->e2 = ed;
ed = EdinTr(tr,p6,p5);
it2->e2 = ed;
if(diag59||diag610) it5->e2 = ed;
}
tr = tr->next;
}
tr = tr4;
for (k = 0; k<4; k++)
{
if (PinTr(p1,tr)) t41 = tr;
else if (PinTr(p2,tr)) t42 = tr;
else if (PinTr(p3,tr)) t43 = tr;
else
{
t44 = tr;
ed = EdinTr(tr,p5,p10);
it1->e3 = ed;
if(diag59) it8->e2 = ed;
else if(diag610) it5->e1 = ed;
ed = EdinTr(tr,p7,p10);
it3->e3 = ed;
if(diag610) it7->e2 = ed;
else if(diag78) it8->e2 = ed;
ed = EdinTr(tr,p5,p7);
it2->e3 = ed;
if(diag59) it7->e2 = ed;
else if(diag78) it6->e1 = ed;
}
tr = tr->next;
}
// four new triangles on the octahedron
SetTP(it1,0);
SetTP(it2,0);
SetTP(it3,0);
SetTP(it4,0);
// four new edges in the octahedron
it5->e3 = ie1;
it6->e3 = ie1;
it7->e3 = ie1;
it8->e3 = ie1;
SetTP(it5,0);
SetTP(it6,0);
SetTP(it7,0);
SetTP(it8,0);
// Build the new tetrahedra.
// tetrahedra generated by chopping off the 4 corners
SetTdP(td1,t21,t31,t41,it1,td);
SetTdE(td1,t21,t31,t41,it1);
neighbourAtFace(td1,n2,t21);
neighbourAtFace(td1,n3,t31);
neighbourAtFace(td1,n4,t41);
SetTdP(td2,t12,t32,t42,it2,td);
SetTdE(td2,t12,t32,t42,it2);
neighbourAtFace(td2,n1,t12);
neighbourAtFace(td2,n3,t32);
neighbourAtFace(td2,n4,t42);
SetTdP(td3,t13,t23,t43,it3,td);
SetTdE(td3,t13,t23,t43,it3);
neighbourAtFace(td3,n1,t13);
neighbourAtFace(td3,n2,t23);
neighbourAtFace(td3,n4,t43);
SetTdP(td4,t14,t24,t34,it4,td);
SetTdE(td4,t14,t24,t34,it4);
neighbourAtFace(td4,n1,t14);
neighbourAtFace(td4,n2,t24);
neighbourAtFace(td4,n3,t34);
// tetrahedra generated by the inner octahedron
if (diag59)
{
SetTdP(td5,t33,it5,it4,it6,td);
SetTdE(td5,t33,it5,it4,it6);
neighbourAtFace(td5,n3,t33);
innerNeighbourAtFace(td5,it5,td6);
innerNeighbourAtFace(td5,it4,td4);
innerNeighbourAtFace(td5,it6,td7);
innerNeighbourAtFace(td4,it4,td5);
SetTdP(td6,t11,it7,it2,it5,td);
SetTdE(td6,t11,it7,it2,it5);
neighbourAtFace(td6,n1,t11);
innerNeighbourAtFace(td6,it7,td8);
innerNeighbourAtFace(td6,it2,td2);
innerNeighbourAtFace(td6,it5,td5);
innerNeighbourAtFace(td2,it2,td6);
SetTdP(td7,t22,it8,it6,it1,td);
SetTdE(td7,t22,it8,it6,it1);
neighbourAtFace(td7,n2,t22);
innerNeighbourAtFace(td7,it8,td8);
innerNeighbourAtFace(td7,it6,td5);
innerNeighbourAtFace(td7,it1,td1);
innerNeighbourAtFace(td1,it1,td7);
SetTdP(td8,t44,it3,it8,it7,td);
SetTdE(td8,t44,it3,it8,it7);
neighbourAtFace(td8,n4,t44);
innerNeighbourAtFace(td8,it3,td3);
innerNeighbourAtFace(td8,it8,td7);
innerNeighbourAtFace(td8,it7,td6);
innerNeighbourAtFace(td3,it3,td8);
}
else if (diag78)
{
SetTdP(td5,t33,it5,it2,it6,td);
SetTdE(td5,t33,it5,it2,it6);
neighbourAtFace(td5,n3,t33);
innerNeighbourAtFace(td5,it5,td6);
innerNeighbourAtFace(td5,it2,td2);
innerNeighbourAtFace(td5,it6,td7);
innerNeighbourAtFace(td2,it2,td5);
SetTdP(td6,t11,it7,it4,it5,td);
SetTdE(td6,t11,it7,it4,it5);
neighbourAtFace(td6,n1,t11);
innerNeighbourAtFace(td6,it7,td8);
innerNeighbourAtFace(td6,it4,td4);
innerNeighbourAtFace(td6,it5,td5);
innerNeighbourAtFace(td4,it4,td6);
SetTdP(td7,t44,it8,it6,it1,td);
SetTdE(td7,t44,it8,it6,it1);
neighbourAtFace(td7,n4,t44);
innerNeighbourAtFace(td7,it8,td8);
innerNeighbourAtFace(td7,it6,td5);
innerNeighbourAtFace(td7,it1,td1);
innerNeighbourAtFace(td1,it1,td7);
SetTdP(td8,t22,it3,it8,it7,td);
SetTdE(td8,t22,it3,it8,it7);
neighbourAtFace(td8,n2,t22);
innerNeighbourAtFace(td8,it3,td3);
innerNeighbourAtFace(td8,it8,td7);
innerNeighbourAtFace(td8,it7,td6);
innerNeighbourAtFace(td3,it3,td8);
}
else
{
SetTdP(td5,t33,it5,it1,it6,td);
SetTdE(td5,t33,it5,it1,it6);
neighbourAtFace(td5,n3,t33);
innerNeighbourAtFace(td5,it5,td7);
innerNeighbourAtFace(td5,it1,td1);
innerNeighbourAtFace(td5,it6,td6);
innerNeighbourAtFace(td1,it1,td5);
SetTdP(td6,t22,it6,it8,it4,td);
SetTdE(td6,t22,it6,it8,it4);
neighbourAtFace(td6,n2,t22);
innerNeighbourAtFace(td6,it6,td5);
innerNeighbourAtFace(td6,it8,td8);
innerNeighbourAtFace(td6,it4,td4);
innerNeighbourAtFace(td4,it4,td6);
SetTdP(td7,t44,it2,it5,it7,td);
SetTdE(td7,t44,it2,it5,it7);
neighbourAtFace(td7,n4,t44);
innerNeighbourAtFace(td7,it2,td2);
innerNeighbourAtFace(td7,it5,td5);
innerNeighbourAtFace(td7,it7,td8);
innerNeighbourAtFace(td2,it2,td7);
SetTdP(td8,t11,it3,it8,it7,td);
SetTdE(td8,t11,it3,it8,it7);
neighbourAtFace(td8,n1,t11);
innerNeighbourAtFace(td8,it3,td3);
innerNeighbourAtFace(td8,it8,td6);
innerNeighbourAtFace(td8,it7,td7);
innerNeighbourAtFace(td3,it3,td8);
}
td1->type = T_WHITE;
td2->type = T_WHITE;
td3->type = T_WHITE;
td4->type = T_WHITE;
td5->type = T_WHITE;
td6->type = T_WHITE;
td7->type = T_WHITE;
td8->type = T_WHITE;
td->firstSon = td1;
td->type = T_RED;
tetList[tetList.h]->add(td1);
tetList[tetList.h]->add(td2);
tetList[tetList.h]->add(td3);
tetList[tetList.h]->add(td4);
tetList[tetList.h]->add(td5);
tetList[tetList.h]->add(td6);
tetList[tetList.h]->add(td7);
tetList[tetList.h]->add(td8);
// destruction of the faces
delete it1;
delete it2;
delete it3;
delete it4;
delete it5;
delete it6;
delete it7;
delete it8;
if ( !(t1->onBoundary()) ) {
delete t11;
delete t12;
delete t13;
delete t14;
}
if ( !(t2->onBoundary()) ) {
delete t21;
delete t22;
delete t23;
delete t24;
}
if ( !(t3->onBoundary()) ) {
delete t31;
delete t32;
delete t33;
delete t34;
}
if ( !(t4->onBoundary()) ) {
delete t41;
delete t42;
delete t43;
delete t44;
}
td->returnTriangles(trs);
return True;
}
//-------------------------------------------------------------------------
int MESH3:: faceNoInTd(TR3 *t, TET3 *tet)
{
PT3 *p1=tet->p1, *p2=tet->p2, *p3=tet->p3, *p4=tet->p4;
if ( PinTr(p2,t) && PinTr(p3,t) && PinTr(p4,t) ) return 1;
if ( PinTr(p1,t) && PinTr(p3,t) && PinTr(p4,t) ) return 2;
if ( PinTr(p1,t) && PinTr(p2,t) && PinTr(p4,t) ) return 3;
if ( PinTr(p1,t) && PinTr(p2,t) && PinTr(p3,t) ) return 4;
return 0;
}
void MESH3:: innerNeighbourAtFace(TET3* td, TR3* t, TET3* neighbourTd)
{
int noN = faceNoInTd(t,td);
if (neighbourTd == 0)
cout << "\n\terror in innerNeighbourAtFace: neighbourTd == 0\n";
switch (noN) {
case 1 :
if (td->n1 != 0)
cout << "\n\terror in innerNeighbourAtFace\n";
td->n1 = neighbourTd; return;
case 2 :
if (td->n2 != 0)
cout << "\n\terror in innerNeighbourAtFace\n";
td->n2 = neighbourTd; return;
case 3 :
if (td->n3 != 0)
cout << "\n\terror in innerNeighbourAtFace\n";
td->n3 = neighbourTd; return;
case 4 :
if (td->n4 != 0)
cout << "\n\terror in innerNeighbourAtFace\n";
td->n4 = neighbourTd; return;
default: break;
}
cout << "\n\n dangerous situation in innerNeighbourAtFace\n";
return;
}
void MESH3:: setNeighbourInTd(TET3* td, PATCH* n, int noOfNeighbour)
{
switch (noOfNeighbour) {
case 1 :
td->n1 = n;
return;
case 2 :
td->n2 = n;
return;
case 3 :
td->n3 = n;
return;
case 4 :
td->n4 = n;
return;
default:
cout << "\n\n error in setNeighbourInTd\n";
return;
}
}
void MESH3:: neighbourAtFace(TET3* td, PATCH* N, TR3* t)
{
TET3 *nb, *nbSon;
int no, no1, no2, noN;
noN = faceNoInTd(t,td);
nb = N->castToTET3();
if ( nb == nil ) // N is boundary face
{
setNeighbourInTd(td,t,noN);
switch (noN) {
case 1 : SetTP(td,t,2,3,4); td->boundP=BOUNDARY; return;
case 2 : SetTP(td,t,1,4,3); td->boundP=BOUNDARY; return;
case 3 : SetTP(td,t,1,2,4); td->boundP=BOUNDARY; return;
case 4 : SetTP(td,t,1,3,2); td->boundP=BOUNDARY; return;
default:
cout << "\n\n error in neighbourAtFace\n";
return;
}
}
nbSon = nb->firstSon; // N is tetrahedron
if (nbSon == nil) {
setNeighbourInTd(td,nb,noN);
if(!nb->mark) {
no = faceNoInTd(t,nb);
if (no) setNeighbourInTd(nb,td,no);
}
return;
}
int iEnd = 0;
if (nb->type == T_RED) iEnd = 8;
else if (nb->type == T_GREEN_1) iEnd = 2;
else if (nb->type == T_GREEN_2A) iEnd = 4;
else if (nb->type == T_GREEN_2B) iEnd = 3;
else if (nb->type == T_GREEN_3) iEnd = 4;
for (int i=1; i<=iEnd; i++)
{
no1 = faceNoInTd(t,nbSon);
if ( no1 )
{
if ( nbSon->firstSon != 0 )
{
int jEnd = 0;
if (nbSon->type == T_RED) jEnd = 8;
else if (nbSon->type == T_GREEN_1) jEnd = 2;
else if (nbSon->type == T_GREEN_2A) jEnd = 4;
else if (nbSon->type == T_GREEN_2B) jEnd = 3;
else if (nbSon->type == T_GREEN_3) jEnd = 4;
TET3* nbSonSon = nbSon->firstSon;
for (int j=1; j<=jEnd; j++)
{
no2 = faceNoInTd(t,nbSonSon);
if ( no2 )
{
setNeighbourInTd(nbSonSon,td,no2);
setNeighbourInTd(td,nbSonSon,noN);
return;
}
else nbSonSon = nbSonSon->next;
}
}
setNeighbourInTd(nbSon,td,no1);
setNeighbourInTd(td,nbSon,noN);
return;
}
else nbSon = nbSon->next;
}
return;
}
//-------------------------------------------------------------------------
int MESH3:: RefineTdBey(TET3 *td)
{
TET3 *td1, *td2, *td3, *td4, *td5, *td6, *td7, *td8;
TR3 *tr=0, *tr1=0, *tr2=0, *tr3=0, *tr4=0,
*it1=0, *it2=0, *it3=0, *it4=0, *it5=0, *it6=0, *it7=0, *it8,
*t11=0, *t12=0, *t13=0, *t14=0,
*t21=0, *t22=0, *t23=0, *t24=0,
*t31=0, *t32=0, *t33=0, *t34=0,
*t41=0, *t42=0, *t43=0, *t44=0,
*t1, *t2, *t3, *t4;
TR3 *trs[1+4];
PT3 *p1 = td->p1,*p2 = td->p2, *p3 = td->p3, *p4 = td->p4,
*p5, *p6, *p7, *p8, *p9, *p10;
EDG3 *e1 = td->e1, *e2 = td->e2, *e3 = td->e3,
*e4 = td->e4, *e5 = td->e5, *e6 = td->e6,
*ie1, *ed;
PATCH *n1=td->n1, *n2=td->n2, *n3=td->n3, *n4=td->n4;
int k;
int depth;
// Regular refinement of a tetrahedron, Method of J. Bey
td->getTriangles(trs);
t1=trs[1]; t2=trs[2]; t3=trs[3]; t4=trs[4];
// Get 8 new tetrahedra, 8 new triangles, 1 new edges
// return False if not enough space
td1 = getTET3(); // new TET3;
td2 = getTET3();
td3 = getTET3();
td4 = getTET3();
td5 = getTET3();
td6 = getTET3();
td7 = getTET3();
td8 = getTET3();
noOfTetrahedra += 7;
it1 = new TR3; // new TR3;
it2 = new TR3;
it3 = new TR3;
it4 = new TR3;
it5 = new TR3;
it6 = new TR3;
it7 = new TR3;
it8 = new TR3;
noOfTriangles += 8;
ie1 = edgAlloc.Get(); // new EDG3;
noOfEdges += 1;
if ((td8==0)||(it8==0)||(ie1==0))
{
cout << " no memory in RefineTd \n";
return False;
}
depth = (td->depth) + 1;
ie1->depth = depth;
if (ie1->depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[ie1->depth]->add(ie1);
it1->depth = depth;
it2->depth = depth;
it3->depth = depth;
it4->depth = depth;
it5->depth = depth;
it6->depth = depth;
it7->depth = depth;
it8->depth = depth;
// Set the "son" triangles tik, i=1,2,3,4; k=1,2,3,4;
// if neeccessary refine the triangles
// Triangle 1
if (!RefineTr(t1))
{
cout << " failed in RefineTr called by RefineTd \n";
return False;
}
setTypeOfInnerEdges(t1);
tr1 = t1->firstSon;
if (tr1==0) return False; // No space for Refinement
// Triangle 2
if (!RefineTr(t2))
{
cout << " failed in RefineTr called by RefineTd \n";
return False;
}
tr2 = t2->firstSon;
setTypeOfInnerEdges(t2);
if (tr2==0) return False; // No space for Refinement
// Triangle 3
if (!RefineTr(t3))
{
cout << " failed in RefineTr called by RefineTd \n";
return False;
}
setTypeOfInnerEdges(t3);
tr3 = t3->firstSon;
if (tr3==0) return False; /* No space for Refinement */
// Triangle 4
if (!RefineTr(t4))
{
cout << " failed in RefineTr called by RefineTd \n";
return False;
}
setTypeOfInnerEdges(t4);
tr4 = t4->firstSon;
if (tr4==0) return False; /* No space for Refinement */
p5 = e3->pm;
p6 = e5->pm;
p7 = e2->pm;
p8 = e6->pm;
p9 = e4->pm;
p10= e1->pm;
// one new edge in the octahedron
ie1->p1 = p6;
ie1->p2 = p10;
ie1->pm = 0;
// new faces
tr = tr1;
for (k = 0; k<4; k++)
{
if (PinTr(p2,tr)) t12 = tr;
else if (PinTr(p3,tr)) t13 = tr;
else if (PinTr(p4,tr)) t14 = tr;
else
{
t11 = tr;
ed = EdinTr(tr,p6,p7);
it2->e1 = ed;
it7->e1 = ed;
ed = EdinTr(tr,p7,p9);
it3->e1 = ed;
ed = EdinTr(tr,p6,p9);
it4->e1 = ed;
it8->e2 = ed;
}
tr = tr->next;
}
tr = tr2;
for (k = 0; k<4; k++)
{
if (PinTr(p1,tr)) t21 = tr;
else if (PinTr(p3,tr)) t23 = tr;
else if (PinTr(p4,tr)) t24 = tr;
else
{
t22 = tr;
ed = EdinTr(tr,p8,p10);
it1->e1 = ed;
it6->e1 = ed;
ed = EdinTr(tr,p8,p9);
it4->e2 = ed;
ed = EdinTr(tr,p9,p10);
it3->e2 = ed;
it8->e1 = ed;
}
tr = tr->next;
}
tr = tr3;
for (k = 0; k<4; k++)
{
if (PinTr(p1,tr)) t31 = tr;
else if (PinTr(p2,tr)) t32 = tr;
else if (PinTr(p4,tr)) t34 = tr;
else
{
t33 = tr;
ed = EdinTr(tr,p5,p8);
it1->e2 = ed;
ed = EdinTr(tr,p6,p8);
it4->e3 = ed;
it6->e2 = ed;
ed = EdinTr(tr,p6,p5);
it2->e2 = ed;
it5->e2 = ed;
}
tr = tr->next;
}
tr = tr4;
for (k = 0; k<4; k++)
{
if (PinTr(p1,tr)) t41 = tr;
else if (PinTr(p2,tr)) t42 = tr;
else if (PinTr(p3,tr)) t43 = tr;
else
{
t44 = tr;
ed = EdinTr(tr,p5,p10);
it1->e3 = ed;
it5->e1 = ed;
ed = EdinTr(tr,p7,p10);
it3->e3 = ed;
it7->e2 = ed;
ed = EdinTr(tr,p5,p7);
it2->e3 = ed;
}
tr = tr->next;
}
// four new triangles on the octahedron
SetTP(it1,0);
SetTP(it2,0);
SetTP(it3,0);
SetTP(it4,0);
// four new triangles in the octahedron
it5->e3 = ie1;
it6->e3 = ie1;
it7->e3 = ie1;
it8->e3 = ie1;
SetTP(it5,0);
SetTP(it6,0);
SetTP(it7,0);
SetTP(it8,0);
// Build the new tetrahedra.
// tetrahedra generated by chopping off the 4 corners
td1->p1 = p1; td1->p2 = p5; td1->p3 = p10; td1->p4 = p8;
SetTdPBey(td1,td);
SetTdE(td1,it1,t21,t31,t41);
neighbourAtFace(td1,n2,t21);
neighbourAtFace(td1,n3,t31);
neighbourAtFace(td1,n4,t41);
td2->p1 = p5; td2->p2 = p2; td2->p3 = p7; td2->p4 = p6;
SetTdPBey(td2,td);
SetTdE(td2,t12,it2,t32,t42);
neighbourAtFace(td2,n1,t12);
neighbourAtFace(td2,n3,t32);
neighbourAtFace(td2,n4,t42);
td3->p1 = p10; td3->p2 = p7; td3->p3 = p3; td3->p4 = p9;
SetTdPBey(td3,td);
SetTdE(td3,t13,t23,it3,t43);
neighbourAtFace(td3,n1,t13);
neighbourAtFace(td3,n2,t23);
neighbourAtFace(td3,n4,t43);
td4->p1 = p8; td4->p2 = p6; td4->p3 = p9; td4->p4 = p4;
SetTdPBey(td4,td);
SetTdE(td4,t14,t24,t34,it4);
neighbourAtFace(td4,n1,t14);
neighbourAtFace(td4,n2,t24);
neighbourAtFace(td4,n3,t34);
// tetrahedra generated by the inner octahedron
td5->p1 = p5; td5->p2 = p10; td5->p3 = p8; td5->p4 = p6;
SetTdPBey(td5,td);
SetTdE(td5,it6,t33,it5,it1);
neighbourAtFace(td5,n3,t33);
td6->p1 = p5; td6->p2 = p10; td6->p3 = p7; td6->p4 = p6;
SetTdPBey(td6,td);
SetTdE(td6,it7,it2,it5,t44);
neighbourAtFace(td6,n4,t44);
td7->p1 = p10; td7->p2 = p8; td7->p3 = p6; td7->p4 = p9;
SetTdPBey(td7,td);
SetTdE(td7,it4,it8,t22,it6);
neighbourAtFace(td7,n2,t22);
td8->p1 = p10; td8->p2 = p7; td8->p3 = p6; td8->p4 = p9;
SetTdPBey(td8,td);
SetTdE(td8,t11,it8,it3,it7);
neighbourAtFace(td8,n1,t11);
innerNeighbourAtFace(td1,it1,td5);
innerNeighbourAtFace(td2,it2,td6);
innerNeighbourAtFace(td3,it3,td8);
innerNeighbourAtFace(td4,it4,td7);
innerNeighbourAtFace(td5,it6,td7);
innerNeighbourAtFace(td5,it5,td6);
innerNeighbourAtFace(td5,it1,td1);
innerNeighbourAtFace(td6,it7,td8);
innerNeighbourAtFace(td6,it2,td2);
innerNeighbourAtFace(td6,it5,td5);
innerNeighbourAtFace(td7,it4,td4);
innerNeighbourAtFace(td7,it8,td8);
innerNeighbourAtFace(td7,it6,td5);
innerNeighbourAtFace(td8,it8,td7);
innerNeighbourAtFace(td8,it3,td3);
innerNeighbourAtFace(td8,it7,td6);
td1->type = T_WHITE;
td2->type = T_WHITE;
td3->type = T_WHITE;
td4->type = T_WHITE;
td5->type = T_WHITE;
td6->type = T_WHITE;
td7->type = T_WHITE;
td8->type = T_WHITE;
td->firstSon = td1;
td->type = T_RED;
tetList[tetList.h]->add(td1);
tetList[tetList.h]->add(td2);
tetList[tetList.h]->add(td3);
tetList[tetList.h]->add(td4);
tetList[tetList.h]->add(td5);
tetList[tetList.h]->add(td6);
tetList[tetList.h]->add(td7);
tetList[tetList.h]->add(td8);
// destruction of the faces
delete it1;
delete it2;
delete it3;
delete it4;
delete it5;
delete it6;
delete it7;
delete it8;
if ( !(t1->onBoundary()) ) {
delete t11;
delete t12;
delete t13;
delete t14;
}
if ( !(t2->onBoundary()) ) {
delete t21;
delete t22;
delete t23;
delete t24;
}
if ( !(t3->onBoundary()) ) {
delete t31;
delete t32;
delete t33;
delete t34;
}
if ( !(t4->onBoundary()) ) {
delete t41;
delete t42;
delete t43;
delete t44;
}
td->returnTriangles(trs);
return True;
}
//-------------------------------------------------------------------------
EDG3* MESH3:: findEdgInNeighbour(PT3* p1, PT3* p2, TET3* td)
{
TET3 *tdSon;
EDG3 *e;
int iEnd;
tdSon = td->firstSon;
if (tdSon == nil) {
cout << "\n\n error in findEdgInNeighbour: tdSon == nil\n" ;
return 0;
}
iEnd = 8;
if (td->type == T_GREEN_2B) iEnd = 3;
for (int i=1; i<=iEnd; i++)
{
e = EdinTet(p1,p2,tdSon);
if (e != 0) return e;
else tdSon = tdSon->next;
}
return 0;
}
//-------------------------------------------------------------------------
void MESH3:: setTypeOfInnerEdges(TR3 *t)
{
TR3 *tSon4;
tSon4 = t->firstSon->next->next->next;
tSon4->e1->type = T_RED;
tSon4->e2->type = T_RED;
tSon4->e3->type = T_RED;
return;
}
int MESH3:: RefineTr(TR3 *t)
{
TR3 *t1, *t2, *t3, *t4;
EDG3 *e1 = t->e1, *e2 = t->e2, *e3 = t->e3, *ereg,
*ie1,*ie2,*ie3, *e11, *e12, *e21, *e22, *e31, *e32;
TET3 *nt;
int depth = (t->depth)+1;
Bool nbRefined, nbFits;
// Get 4 new triangles, 3 new edges, return False if not enough space
t1 = new TR3;
t2 = new TR3;
t3 = new TR3;
t4 = new TR3;
// Set the "son" edges eik, i=1,2,3; k=1,2;
// if neeccessary refine the edges
if ((e1->firstSon)==0) /* Edge not refined */
{
e11 = RefineEdg(e1);
if (e11==0) return False; /* No space for Refinement */
}
else e11 = e1->firstSon; /* Take the existing "son" edge */
e12 = e11->next;
if ((e11->p1)==(t->p3))
{ ereg = e11; e11 = e12; e12 = ereg; }
if ((e2->firstSon)==0)
{
e21 = RefineEdg(e2);
if (e21==0) return False;
}
else e21 = e2->firstSon;
e22 = e21->next;
if ((e21->p1)==(t->p1))
{ ereg = e21; e21 = e22; e22 = ereg; }
if ((e3->firstSon)==0)
{
e31 = RefineEdg(e3);
if (e31==0) return False;
}
else e31 = e3->firstSon;
e32 = e31->next;
if ((e31->p1)==(t->p2))
{ ereg = e31; e31 = e32; e32 = ereg; }
ie1 = edgAlloc.Get();
ie2 = edgAlloc.Get();
ie3 = edgAlloc.Get();
if ((t4==0)||(ie3==0)) return False;
ie1->depth = ie2->depth = ie3->depth = depth;
ie1->boundP = t->boundP;
ie2->boundP = t->boundP;
ie3->boundP = t->boundP;
ie1->classA = t->classA;
ie2->classA = t->classA;
ie3->classA = t->classA;
ie1->isoP = t->isoP;
ie2->isoP = t->isoP;
ie3->isoP = t->isoP;
// Build the new triangles.
InnerEdges(t,e1,e2,e3,ie1,ie2,ie3);
nt = t->t32; // looking for the second neighbour
nbFits = False; nbRefined = False;
if ( nt != 0 )
{
if ( PinTet(t->p1,nt) && PinTet(t->p2,nt) && PinTet(t->p3,nt) )
nbFits = True;
if ( nt->refined() ) nbRefined=True;
}
if ( nbFits && nbRefined )
{
EDG3 *ne1, *ne2, *ne3;
ne1 = findEdgInNeighbour(ie1->p1,ie1->p2,nt);
ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
ne3 = findEdgInNeighbour(ie3->p1,ie3->p2,nt);
edgAlloc.Return(ie1);
edgAlloc.Return(ie2);
edgAlloc.Return(ie3);
ie1=ne1; ie2=ne2; ie3=ne3;
}
else
{
noOfEdges += 3;
noOfTriangles += 3;
if (depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[depth]->add(ie1);
edgList[depth]->add(ie2);
edgList[depth]->add(ie3);
}
t1->e1 = e32; t1->e2 = e11; t1->e3 = ie1;
t2->e1 = e12; t2->e2 = e21; t2->e3 = ie2;
t3->e1 = e22; t3->e2 = e31; t3->e3 = ie3;
t4->e3 = ie1; t4->e2 = ie3; t4->e1 = ie2;
SetTP(t1,t);
SetTP(t2,t);
SetTP(t3,t);
SetTP(t4,t);
t->firstSon = t1;
t->type = T_RED;
t1->next = t2; t2->next = t3, t3->next = t4;
t2->prev = t1; t3->prev = t2, t4->prev = t3;
if (t->boundP != INTERIOR)
{
if (depth > trList.h) // create new list at depth if necessary
{
DList<TR3>* trL = new DList<TR3>;
trList.push(trL);
}
trList[depth]->add(t1);
trList[depth]->add(t2);
trList[depth]->add(t3);
trList[depth]->add(t4);
// noOfTriangles += 3;
}
return True;
}
//-------------------------------------------------------------------------
EDG3* MESH3:: RefineEdg(EDG3 *ed)
{
EDG3 *e1, *e2;
PT3 *pm = 0;
int depth;
// Refinement of an edge.
// -- Get two new edges; return False if not enough space
e1 = edgAlloc.Get();
e2 = edgAlloc.Get();
if (e2==0) return 0;
noOfEdges++;
pm = ptAlloc.Get();
pm->reset();
noOfPoints++;
// Call the user midpoint routine and set the new fields
pm->classA = ed->classA;
pm->boundP = ed->boundP;
pm->depth = ed->depth + 1;
if ( (ed->isoP == CIRCULAR) && (ed->boundP != INTERIOR) )
setXYZpmCircle(pm, ed);
else
if ( (ed->boundP != INTERIOR) && (ed->isoP == CYLINDRICAL) )
setXYZpmCyl(pm, ed);
else
{
pm->x = HALF*(ed->p1->x + ed->p2->x);
pm->y = HALF*(ed->p1->y + ed->p2->y);
pm->z = HALF*(ed->p1->z + ed->p2->z);
}
e1->p1 = ed->p1; e1->p2 = pm;
e1->boundP = ed->boundP; e1->father = ed; e1->type = ed->type;
e2->p1 = pm; e2->p2 = ed->p2;
e2->boundP = ed->boundP; e2->father = ed; e2->type = ed->type;
e1->classA = e2->classA = ed->classA;
e1->isoP = e2->isoP = ed->isoP;
ed->firstSon = e1; ed->pm =pm;
e1->depth = e2->depth = ed->depth+1;
depth = e1->depth;
if (depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[depth]->add(e1);
edgList[depth]->add(e2);
if (pm->depth > ptList.h)
{
DList<PT3>* ptL = new DList<PT3>;
ptList.push(ptL);
}
ptList[pm->depth]->add(pm);
return e1;
}
//-------------------------------------------------------------------------
void MESH3:: SetTP(TR3 *t,TR3 *f)
{
PT3 *P11 = (t->e1)->p1, *P12 = (t->e1)->p2,
*P21 = (t->e2)->p1, *P22 = (t->e2)->p2,
*P31 = (t->e3)->p1, *P32 = (t->e3)->p2;
// set the points of a triangle, computed from the edges
t->p1 = ((P31==P11)||(P31==P12))?P32:P31;
t->p2 = ((P11==P21)||(P11==P22))?P12:P11;
t->p3=((P21==P31)||(P21==P32))?P22:P21;
t->father = f;
if (f!=0)
{
t->depth = (f->depth)+1;
t->classA = f->classA;
t->isoP = f->isoP;
t->boundP = f->boundP;
}
else
t->boundP = INTERIOR;
if ( (t->depth) > maxDepth )
maxDepth = t->depth;
return;
}
//-------------------------------------------------------------------------
void MESH3:: SetTP(TET3 *tet, TR3 *t, int noP1, int noP2, int noP3)
{
PT3 *tetPts[1+4];
tet->getPoints(tetPts);
EDG3 *e1 = t->e1, *e2 = t->e2, *e3 = t->e3;
PT3 *P11 = e1->p1, *P12 = e1->p2,
*P21 = e2->p1, *P22 = e2->p2,
*P31 = e3->p1, *P32 = e3->p2;
t->p1 = tetPts[noP1];
t->p2 = tetPts[noP2];
t->p3 = tetPts[noP3];
if( (P11!=t->p1) && (P12!=t->p1) ) t->e1 = e1;
else if( (P21!=t->p1) && (P22!=t->p1) ) t->e1 = e2;
else if( (P31!=t->p1) && (P32!=t->p1) ) t->e1 = e3;
if( (P11!=t->p2) && (P12!=t->p2) ) t->e2 = e1;
else if( (P21!=t->p2) && (P22!=t->p2) ) t->e2 = e2;
else if( (P31!=t->p2) && (P32!=t->p2) ) t->e2 = e3;
if( (P11!=t->p3) && (P12!=t->p3) ) t->e3 = e1;
else if( (P21!=t->p3) && (P22!=t->p3) ) t->e3 = e2;
else if( (P31!=t->p3) && (P32!=t->p3) ) t->e3 = e3;
return;
}
//-------------------------------------------------------------------------
void TET3:: completeT(TR3 *t,TR3 *ft) const
{
/*
PT3 *P11 = (t->e1)->p1, *P12 = (t->e1)->p2,
*P21 = (t->e2)->p1, *P22 = (t->e2)->p2,
*P31 = (t->e3)->p1, *P32 = (t->e3)->p2;
t->p1 = ((P31==P11)||(P31==P12))?P32:P31;
t->p2 = ((P11==P21)||(P11==P22))?P12:P11;
t->p3=((P21==P31)||(P21==P32))?P22:P21;
*/
t->father = ft;
if (ft!=0)
{
t->depth = (ft->depth)+1;
t->classA = ft->classA;
t->isoP = ft->isoP;
t->boundP = ft->boundP;
}
else
t->boundP = INTERIOR;
return;
}
//-------------------------------------------------------------------------
int MESH3:: TdVolume(TET3 *tet)
{
Real x1, x2, x3, y1, y2, y3, z1, z2, z3;
Real det;
PT3 *P1, *P2, *P3, *P4;
// computes the Jacobian and the orientation of the edges
P1 = tet->p1;
P2 = tet->p2;
P3 = tet->p3;
P4 = tet->p4;
x1 = (P2->x) - (P1->x);
x2 = (P2->y) - (P1->y);
x3 = (P2->z) - (P1->z);
y1 = (P3->x) - (P1->x);
y2 = (P3->y) - (P1->y);
y3 = (P3->z) - (P1->z);
z1 = (P4->x) - (P1->x);
z2 = (P4->y) - (P1->y);
z3 = (P4->z) - (P1->z);
det = x1*(y2*z3 - y3*z2) - x2*(y1*z3 -y3*z1) + x3*(y1*z2 - y2*z1);
tet->detJ = det; // the volume of the td: fabs(det)*SIXTH ;
return tet->orient();
}
//-------------------------------------------------------------------------
void MESH3:: SetTdP(TET3 *td, TR3* t1, TR3* t2, TR3* t3, TR3* t4, TET3 *f)
{
PT3 *P11 = t1->p1, *P12 = t1->p2, *P13 = t1->p3,
*P21 = t2->p1, *P22 = t2->p2, *P23 = t2->p3,
*P31 = t3->p1, *P32 = t3->p2, *P33 = t3->p3,
*P41 = t4->p1, *P42 = t4->p2, *P43 = t4->p3;
PT3 *P2, *P3;
// set the points of a tetrahedron td, computed from the faces;
// td inherits properties of its father f
if ((P31!=P11)&&(P31!=P12)&&(P31!=P13))
td->p1 = P31;
else if ((P32!=P11)&&(P32!=P12)&&(P32!=P13))
td->p1 = P32;
else td->p1 = P33;
if ((P41!=P21)&&(P41!=P22)&&(P41!=P23))
td->p2 = P41;
else if ((P42!=P21)&&(P42!=P22)&&(P42!=P23))
td->p2 = P42;
else td->p2 = P43;
if ((P11!=P31)&&(P11!=P32)&&(P11!=P33))
td->p3 = P11;
else if ((P12!=P31)&&(P12!=P32)&&(P12!=P33))
td->p3 = P12;
else td->p3 = P13;
if ((P21!=P41)&&(P21!=P42)&&(P21!=P43))
td->p4 = P21;
else if ((P22!=P41)&&(P22!=P42)&&(P22!=P43))
td->p4 = P22;
else td->p4 = P23;
TdVolume(td);
if (rightHand)
{
P2 = td->p2;
P3 = td->p3;
if (!(td->orient()))
{
td->p2 = P3;
td->p3 = P2;
td->detJ = Abs(td->detJ); // td->orient = True;
}
}
td->n1 = nil; td->n2 = nil; td->n3 = nil; td->n4 = nil;
td->father = f;
if (f!=0)
{
td->depth = f->depth +1;
td->classA = f->classA;
}
if (td->depth > maxDepth) maxDepth = td->depth;
return;
}
//-------------------------------------------------------------------------
void MESH3:: SetTdPBey(TET3 *td,TET3 *f)
{
// set the points of a tetrahedron td, computed from the faces;
// td inherits properties of its father f
TdVolume(td);
td->n1 = nil; td->n2 = nil; td->n3 = nil; td->n4 = nil;
td->father = f;
if (f!=0)
{
td->depth = (f->depth)+1;
td->classA = f->classA;
}
if ( (td->depth)> maxDepth )
maxDepth = td->depth;
return;
}
//-------------------------------------------------------------------------
void MESH3:: SetTPts(TET3 *td, TET3 *f)
{
PT3 *P2, *P3;
// computes volume and orientation in the tetrahedron td;
// td inherits properties of its father f
TdVolume(td); // determinant
if (rightHand)
{
P2 = td->p2;
P3 = td->p3;
if (!(td->orient()))
{
td->p2 = P3;
td->p3 = P2;
td->detJ = Abs(td->detJ); // td->orient() = True;
}
}
td->father = f;
if (f != 0)
{
td->depth = (f->depth)+1;
td->classA = f->classA;
}
if ( (td->depth)> maxDepth )
maxDepth = td->depth;
return;
}
//-------------------------------------------------------------------------
EDG3* MESH3:: FindTdE(TET3 */*td*/,TR3 *t1,TR3 *t2,TR3 *t3,TR3 *t4,PT3 *p1,PT3 *p2)
{
EDG3 *E11, *E12, *E13,
*E21, *E22, *E23,
*E31, *E32, *E33,
*E41, *E42, *E43;
E11 = t1->e1;
if ( ((E11->p1==p1)&&(E11->p2==p2))||((E11->p1==p2)&&(E11->p2==p1)))
return E11;
E12 = t1->e2;
if ( ((E12->p1==p1)&&(E12->p2==p2))||((E12->p1==p2)&&(E12->p2==p1)))
return E12;
E13 = t1->e3;
if ( ((E13->p1==p1)&&(E13->p2==p2))||((E13->p1==p2)&&(E13->p2==p1)))
return E13;
E21 = t2->e1;
if ( ((E21->p1==p1)&&(E21->p2==p2))||((E21->p1==p2)&&(E21->p2==p1)))
return E21;
E22 = t2->e2;
if ( ((E22->p1==p1)&&(E22->p2==p2))||((E22->p1==p2)&&(E22->p2==p1)))
return E22;
E23 = t2->e3;
if ( ((E23->p1==p1)&&(E23->p2==p2))||((E23->p1==p2)&&(E23->p2==p1)))
return E23;
E31 = t3->e1;
if ( ((E31->p1==p1)&&(E31->p2==p2))||((E31->p1==p2)&&(E31->p2==p1)))
return E31;
E32 = t3->e2;
if ( ((E32->p1==p1)&&(E32->p2==p2))||((E32->p1==p2)&&(E32->p2==p1)))
return E32;
E33 = t3->e3;
if ( ((E33->p1==p1)&&(E33->p2==p2))||((E33->p1==p2)&&(E33->p2==p1)))
return E33;
E41 = t4->e1;
if ( ((E41->p1==p1)&&(E41->p2==p2))||((E41->p1==p2)&&(E41->p2==p1)))
return E41;
E42 = t4->e2;
if ( ((E42->p1==p1)&&(E42->p2==p2))||((E42->p1==p2)&&(E42->p2==p1)))
return E42;
E43 = t4->e3;
if ( ((E43->p1==p1)&&(E43->p2==p2))||((E43->p1==p2)&&(E43->p2==p1)))
return E43;
return 0;
}
//-------------------------------------------------------------------------
void MESH3:: SetTdE(TET3 *td,TR3 *t1,TR3 *t2,TR3 *t3,TR3 *t4)
{
PT3 *P1 = td->p1, *P2 = td->p2, *P3 = td->p3, *P4 = td->p4;
// set the edges of a tetrahedron td
td->e1 = FindTdE(td,t1,t2,t3,t4,P1,P3);
td->e2 = FindTdE(td,t1,t2,t3,t4,P2,P3);
td->e3 = FindTdE(td,t1,t2,t3,t4,P1,P2);
td->e4 = FindTdE(td,t1,t2,t3,t4,P3,P4);
td->e5 = FindTdE(td,t1,t2,t3,t4,P2,P4);
td->e6 = FindTdE(td,t1,t2,t3,t4,P1,P4);
return;
}
//-------------------------------------------------------------------------
void MESH3:: InnerEdges(TR3 *t,EDG3 *e1,EDG3 *e2,EDG3 *e3,EDG3 *ie1,
EDG3 *ie2,EDG3 *ie3)
{
// sets the points of the inner edges *ie1,*ie2,*ie3;
// Edges of the triangle t to be refined: *e1,*e2,*e3
ie1->p1 = e1->pm;
ie1->p2 = e3->pm;
ie1->type = T_GREEN_3;
ie1->classA = t->classA; // !!!
ie1->isoP = t->isoP;
ie1->boundP = t->boundP;
ie2->p1 = e2->pm;
ie2->p2 = e1->pm;
ie2->type = T_GREEN_3;
ie2->classA = t->classA;
ie2->isoP = t->isoP;
ie2->boundP = t->boundP;
ie3->p1 = e3->pm;
ie3->p2 = e2->pm;
ie3->type = T_GREEN_3;
ie3->classA = t->classA;
ie3->isoP = t->isoP;
ie3->boundP = t->boundP;
return;
}
//-------------------------------------------------------------------------
int MESH3:: MkGreen(TR3 *t)
{
PT3 *p1, *p2;
EDG3 *e1 = t->e1,*e2 = t->e2,*e3 = t->e3, *eNew=0;
TR3 *t1New=0, *t2New=0;
TET3 *nt;
int depth;
Bool nbRefined, nbFits;
// generates two green triangles as sons of t
if (t->firstSon != 0) return True;
// Testing if one neighbor triangle is refined and setting pi, ei
// corresponding to this neighbor
if ( (e1->type == T_GREEN_1) || (e2->type == T_GREEN_1) ||
(e3->type == T_GREEN_1) )
cout << "Alarm in MkGreen: one edge already green \n";
if ((e1->firstSon)!=0)
{
p1 = t->p1; p2 = t->p2;
e1 = t->e1; e2 = t->e2; e3 = t->e3;
}
else
{
if ((e2->firstSon)!=0)
{
p1 = t->p2; p2 = t->p3;
e1 = t->e2; e2 = t->e3; e3 = t->e1;
}
else
{
if ((e3->firstSon)!=0)
{
p1 = t->p3; p2 = t->p1;
e1 = t->e3; e2 = t->e1; e3 = t->e2;
}
else return True; /* nothing to refine */
}
}
// Get one new edge and two new triangles
eNew = edgAlloc.Get();
t1New = new TR3;
t2New = new TR3;
// Set the fields for the new objects
eNew->p1 = p1; eNew->p2 = (e1->firstSon)->p2;
nt = t->t32; // looking for the second neighbour
nbFits = False; nbRefined = False;
if ( nt != 0 )
{
if ( PinTet(t->p1,nt) && PinTet(t->p2,nt) && PinTet(t->p3,nt) )
nbFits = True;
if ( nt->refined() ) nbRefined=True;
}
if ( nbFits && nbRefined )
{
EDG3 *ne1;
ne1 = findEdgInNeighbour(eNew->p1,eNew->p2,nt);
edgAlloc.Return(eNew);
eNew=ne1;
}
else
{
noOfEdges++;
noOfTriangles++;
eNew->type = T_GREEN_1;
eNew->classA = t->classA;
eNew->isoP = t->isoP;
eNew->boundP = t->boundP;
depth = t->depth+1;
eNew->depth = depth;
if (depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[depth]->add(eNew);
}
t1New->e1 = e3; t1New->e3 = eNew; t1New->type = T_GREEN_1;
t2New->e1 = e2; t2New->e2 = eNew; t2New->type = T_GREEN_1;
if ((e1->p1)==p2) // test for orientation
{
t1New->e2 = e1->firstSon; t2New->e3 = (e1->firstSon)->next;
}
else
{
t1New->e2 = (e1->firstSon)->next; t2New->e3 = e1->firstSon;
}
SetTP(t1New,t); SetTP(t2New,t);
t->firstSon = t1New;
t->type = T_GREEN_1;
// Append two new ones at their depth list
depth = t1New->depth;
t1New->next = t2New; t2New->prev = t1New;
if (t->boundP != INTERIOR)
{
if ( depth > trList.h ) // create new list at depth if necessary
{
DList<TR3>* trL = new DList<TR3>;
trList.push(trL);
}
trList[depth]->add(t1New);
trList[depth]->add(t2New);
// noOfTriangles++;
}
return True;
}
//-------------------------------------------------------------------------
int MESH3:: MkGreenIrr(TR3 *t)
{
TR3 *t1=0, *t2=0, *t3=0;
EDG3 *e1 = t->e1, *e2 = t->e2, *e3 = t->e3, *ereg,
*ie1=0,*ie2=0, *e11=0, *e12=0, *e21=0, *e22=0, *e31=0, *e32=0;
TET3 *nt;
int depth;
Bool nbFits, nbRefined, edgeChanged=False;
//cout << "\n MkGreenIrr \n";
// generates three green triangles as sons of t
if (t->firstSon != 0) return True;
// Get 3 new triangles, 2 new edges
t1 = new TR3;
t2 = new TR3;
t3 = new TR3;
ie1 = edgAlloc.Get();
ie2 = edgAlloc.Get();
ie1->depth = t->depth+1;
ie2->depth = t->depth+1;
ie1->next = ie2;
ie2->prev = ie1;
// boundary value on the new triangles and edges
ie1->boundP = t->boundP;
ie1->type = T_GREEN_2B;
ie1->classA = t->classA;
ie1->isoP = t->isoP;
ie2->boundP = t->boundP;
ie2->type = T_GREEN_2B;
ie2->classA = t->classA;
ie2->isoP = t->isoP;
// Set the "son" edges eik, i=1,2,3; k=1,2;
// if neeccessary refine the edges
nt = t->t32; // looking for the second neighbour
nbFits = False; nbRefined = False;
if ( nt != 0 )
{
if ( PinTet(t->p1,nt) && PinTet(t->p2,nt) && PinTet(t->p3,nt) )
nbFits = True;
if ( nt->refined() ) nbRefined=True;
}
if ((e1->firstSon)==0) /* Edge e1 not refined */
{
ie1->p1 = e3->pm;
ie1->p2 = e2->pm;
ie2->p1 = e3->pm;
ie2->p2 = t->p3;
if ( nbFits && nbRefined )
{
EDG3 *ne1, *ne2;
ne1 = findEdgInNeighbour(ie1->p1,ie1->p2,nt);
ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
if ( ne2 == 0 )
{
ie2->p1 = e2->pm;
ie2->p2 = t->p2;
ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
edgeChanged = True;
}
edgAlloc.Return(ie1); edgAlloc.Return(ie2);
ie1=ne1; ie2=ne2;
}
else
{
noOfEdges += 2;
noOfTriangles += 2;
depth = ie1->depth;
if (depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[depth]->add(ie1);
edgList[depth]->add(ie2);
}
e21 = e2->firstSon;
e22 = e21->next;
if ( (e21->p1) == (t->p3) )
{ereg=e21; e21 = e22; e22 = ereg;}
e31 = e3->firstSon;
e32 = e31->next;
if ( (e31->p1) == (t->p2) )
{ereg=e31; e31 = e32; e32 = ereg;}
if (edgeChanged) { EDG3* e = e22; e22 = e32; e32 = e;}
t1->e1 = e31; t1->e2 = ie1; t1->e3 = e21;
t2->e1 = e1; t2->e2 = ie2; t2->e3 = e32;
t3->e1 = ie1; t3->e2 = ie2; t3->e3 = e22;
}
else
if ((e2->firstSon)==0) /* Edge e2 not refined */
{
ie1->p1 = e1->pm;
ie1->p2 = e3->pm;
ie2->p1 = e1->pm;
ie2->p2 = t->p1;
if ( nbFits && nbRefined )
{
EDG3 *ne1, *ne2;
ne1 = findEdgInNeighbour(ie1->p1,ie1->p2,nt);
ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
if ( ne2 == 0 )
{
ie2->p1 = e3->pm;
ie2->p2 = t->p3;
ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
edgeChanged = True;
}
edgAlloc.Return(ie1); edgAlloc.Return(ie2);
ie1=ne1; ie2=ne2;
}
else
{
noOfEdges += 2;
noOfTriangles += 2;
depth = ie1->depth;
if (depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[depth]->add(ie1);
edgList[depth]->add(ie2);
}
e11 = e1->firstSon;
e12 = e11->next;
if ( (e11->p1) == (t->p3) )
{ereg=e11; e11 = e12; e12 = ereg;}
e31 = e3->firstSon;
e32 = e31->next;
if ( (e31->p1) == (t->p2) )
{ereg=e31; e31 = e32; e32 = ereg;}
if (edgeChanged) { EDG3* e = e12; e12 = e31; e31 = e;}
t1->e1 = e32; t1->e2 = e11; t1->e3 = ie1;
t2->e1 = e12; t2->e2 = e2; t2->e3 = ie2;
t3->e1 = e31; t3->e2 = ie1; t3->e3 = ie2;
}
else
if ((e3->firstSon)==0) /* Edge e3 not refined */
{
ie1->p1 = e2->pm;
ie1->p2 = e1->pm;
ie2->p1 = e2->pm;
ie2->p2 = t->p2;
if ( nbFits && nbRefined )
{
EDG3 *ne1, *ne2;
ne1 = findEdgInNeighbour(ie1->p1,ie1->p2,nt);
ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
if ( ne2 == 0 )
{
ie2->p1 = e1->pm;
ie2->p2 = t->p1;
ne2 = findEdgInNeighbour(ie2->p1,ie2->p2,nt);
edgeChanged = True;
}
edgAlloc.Return(ie1); edgAlloc.Return(ie2);
ie1=ne1; ie2=ne2;
}
else
{
noOfEdges += 2;
noOfTriangles += 2;
depth = ie1->depth;
if (depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[depth]->add(ie1);
edgList[depth]->add(ie2);
}
e11 = e1->firstSon;
e12 = e11->next;
if ( (e11->p1) == (t->p3) )
{ereg=e11; e11 = e12; e12 = ereg;}
e21 = e2->firstSon;
e22 = e21->next;
if ( (e21->p1) == (t->p3) )
{ereg=e21; e21 = e22; e22 = ereg;}
if (edgeChanged) { EDG3* e = e11; e11 = e21; e21 = e;}
t1->e1 = e12; t1->e2 = e22; t1->e3 = ie1;
t2->e1 = ie2; t2->e2 = e21; t2->e3 = e3;
t3->e1 = e11; t3->e2 = ie1; t3->e3 = ie2;
}
// Build the new triangles.
t1->type = T_GREEN_2B; SetTP(t1,t);
t2->type = T_GREEN_2B; SetTP(t2,t);
t3->type = T_GREEN_2B; SetTP(t3,t);
t->firstSon = t1;
t->type = T_GREEN_2B;
// Append two new ones at their depth list
depth = t1->depth;
t1->next = t2; t2->next = t3;
t2->prev = t1; t3->prev = t2;
if (t->boundP != INTERIOR)
{
if (depth > trList.h) // create new list at depth if necessary
{
DList<TR3>* trL = new DList<TR3>;
trList.push(trL);
}
trList[depth]->add(t1);
trList[depth]->add(t2);
trList[depth]->add(t3);
// noOfTriangles += 2;
}
return True;
}
//------------------------------------------------------------------------
int MESH3:: MarkGreenTd(TET3 *td)
{
EDG3 *eds[1+6], *ed;
TR3 *trs[1+4];
int i, noOfRefEdg = 0;
// marks green tetrahedra to be refined
if ( (td->type) != T_WHITE )
cout << " Alarm in MarkGreenTd: td not white \n";
td->getEdges(eds);
for (i=1; i<=6; i++)
{
ed = eds[i];
if ( (ed->firstSon) != 0 ) noOfRefEdg++;
}
if ( noOfRefEdg == 1)
{
td->mark = T_GREEN_1;
moreGreen1Td++;
}
else
if ( noOfRefEdg == 2)
{
td->mark = T_GREEN_2A;
moreGreen2ATd++;
td->getTriangles(trs);
for ( i=1; i <= 4; i++)
{
if ( TestFace(trs[i]) == 2 )
{
td->mark = T_GREEN_2B;
moreGreen2ATd--;
moreGreen2BTd++;
}
}
td->returnTriangles(trs);
}
else
if ( noOfRefEdg == 3)
{
td->mark = T_GREEN_3;
moreGreen3Td++;
}
else td->mark = False;
return True;
}
//-------------------------------------------------------------------------
int MESH3:: MkGreen1Td(TET3 *td)
{
PT3 *p1=0, *p2=0, *p3=0, *p4=0, *pm=0;
TR3 *t=0, *t1=0, *t2=0, *t3=0, *t4=0;
TR3 *t1New=0, *t21=0, *t22=0, *t41=0, *t42=0;
TR3 *trs[1+4];
TET3 *td1,*td2;
PATCH *n1=0, *n2=0, *n3=0, *n4=0;
EDG3 *e1=td->e1,*e2=td->e2,*e3=td->e3,*e4=td->e4, *e5=td->e5, *e6=td->e6;
EDG3 *ed, *e1New=0, *e2New=0, *e3New=0;
int refEdges=0, depth = (td->depth)+1;
// refinement of td into two green tetrahedra
if (td->mark == T_GREEN_1) td->mark = False;
else return True;
td->getTriangles(trs);
// Get one new triangle and two tetrahedra,
// if not enough space: return False
t1New = new TR3;
noOfTriangles++;
t1New->depth = depth;
td1 = getTET3();
td2 = getTET3();
noOfTetrahedra++;
// Testing which edge is refined
ed = e1;
if ( (ed->firstSon)!=0 )
{
refEdges +=1;
p1=td->p1; p2=td->p2; p3=td->p3; p4=td->p4;
pm = ed->pm; e3New = e5;
t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
n1 = td->n1; n2 = td->n2; n3 = td->n3; n4 = td->n4;
}
ed = e2;
if ( (ed->firstSon)!=0 )
{
refEdges +=1;
if (refEdges > 1 )
{
cout << "\n\n alarm in MkGreenI \n";
return False;
}
p1=td->p3; p2=td->p1; p3=td->p2; p4=td->p4;
pm = ed->pm; e3New = e6;
t1 = trs[3]; t2 = trs[1]; t3 = trs[2]; t4 = trs[4];
n1 = td->n3; n2 = td->n1; n3 = td->n2; n4 = td->n4;
}
ed = e3;
if ( (ed->firstSon)!=0 )
{
refEdges +=1;
if (refEdges > 1 )
{
cout << "\n\n alarm in MkGreenI \n";
return False;
}
p1=td->p2; p2=td->p3; p3=td->p1; p4=td->p4;
pm = ed->pm; e3New = e4;
t1 = trs[2]; t2 = trs[3]; t3 = trs[1]; t4 = trs[4];
n1 = td->n2; n2 = td->n3; n3 = td->n1; n4 = td->n4;
}
ed = e4;
if ( (ed->firstSon)!=0 )
{
refEdges +=1;
if (refEdges > 1 )
{
cout << "\n\n alarm in MkGreenI \n";
return False;
}
p1=td->p4; p2=td->p1; p3=td->p3; p4=td->p2;
pm = ed->pm; e3New = e3;
t1 = trs[4]; t2 = trs[1]; t3 = trs[3]; t4 = trs[2];
n1 = td->n4; n2 = td->n1; n3 = td->n3; n4 = td->n2;
}
ed = e5;
if ( (ed->firstSon)!=0 )
{
refEdges +=1;
if (refEdges > 1 )
{
cout << "\n\n alarm in MkGreenI \n";
return False;
}
p1=td->p2; p2=td->p3; p3=td->p4; p4=td->p1;
pm = ed->pm; e3New = e1;
t1 = trs[2]; t2 = trs[3]; t3 = trs[4]; t4 = trs[1];
n1 = td->n2; n2 = td->n3; n3 = td->n4; n4 = td->n1;
}
ed = e6;
if ( (ed->firstSon)!=0 )
{
refEdges +=1;
if (refEdges > 1 )
{
cout << "\n\n alarm in MkGreenI \n";
return False;
}
p1=td->p1; p2=td->p3; p3=td->p4; p4=td->p2;
pm = ed->pm; e3New = e2;
t1 = trs[1]; t2 = trs[3]; t3 = trs[4]; t4 = trs[2];
n1 = td->n1; n2 = td->n3; n3 = td->n4; n4 = td->n2;
}
// Set the fields for the new objects
MkGreen(t2);
t = t2->firstSon;
if ( PinTr(p1,t) )
{ t21 = t; t22 = t->next; }
else
{ t22 = t; t21 = t->next; }
e1New = EdinTr(t,p4,pm);
MkGreen(t4);
t = t4->firstSon;
if ( PinTr(p1,t) )
{ t41 = t; t42 = t->next; }
else
{ t42 = t; t41 = t->next; }
e2New = EdinTr(t,p2,pm);
t1New->e1 = e1New; t1New->e2 = e2New; t1New->e3 = e3New;
SetTP(t1New,0); t1New->type = T_GREEN_3;
SetTdP(td1,t21,t3,t41,t1New,td);
SetTdE(td1,t21,t3,t41,t1New);
neighbourAtFace(td1,n2,t21);
neighbourAtFace(td1,n3,t3);
neighbourAtFace(td1,n4,t41);
td1->type = T_GREEN_1;
SetTdP(td2,t22,t1,t42,t1New,td);
SetTdE(td2,t22,t1,t42,t1New);
neighbourAtFace(td2,n2,t22);
neighbourAtFace(td2,n1,t1);
neighbourAtFace(td2,n4,t42);
td2->type = T_GREEN_1;
innerNeighbourAtFace(td1,t1New,td2);
innerNeighbourAtFace(td2,t1New,td1);
td->firstSon = td1;
td->type = T_GREEN_1;
tetList[tetList.h]->add(td1);
tetList[tetList.h]->add(td2);
// destruction of the faces
delete t1New;
if ( !(t2->onBoundary()) ) {
delete t21;
delete t22;
}
if ( !(t4->onBoundary()) ) {
delete t41;
delete t42;
}
td->returnTriangles(trs);
return True;
}
//-------------------------------------------------------------------------
int MESH3:: MkGreen2ATd(TET3 *td)
{
PT3 *p1=0, *p2=0, *p3=0, *p4=0, *pm1=0, *pm2=0;
TR3 *t=0, *t1=0, *t2=0, *t3=0, *t4=0;
TR3 *t1New=0, *t2New=0, *t3New=0, *t4New=0;
TR3 *t11=0, *t12=0, *t21=0, *t22=0, *t31=0, *t32=0, *t41=0, *t42=0;
TR3 *trs[1+4];
TET3 *td1=0, *td2=0, *td3=0, *td4=0;
PATCH *n1=0, *n2=0, *n3=0, *n4=0;
EDG3 *e1=td->e1,*e2=td->e2,*e3=td->e3,*e4=td->e4, *e5=td->e5, *e6=td->e6;
EDG3 *ed, *edNew=0;
EDG3 *e1New, *e2New, *e3New, *e4New, *e5New, *e6New, *e7New, *e8New;
int refEdges, depth = (td->depth)+1;
// Refinement of td into four green tetrahedra
if ( (td->mark == T_GREEN_2A) ) td->mark = False;
else return True;
td->getTriangles(trs);
// Get one new edge, four new triangles and four tetrahedra,
// if not enough space: return False
edNew = edgAlloc.Get();
noOfEdges += 1;
edNew->depth = depth;
if (depth > edgList.h) // create new list at depth if necessary
{
DList<EDG3>* edgL = new DList<EDG3>;
edgList.push(edgL);
}
edgList[depth]->add(edNew);
t1New = new TR3;
t2New = new TR3;
t3New = new TR3;
t4New = new TR3;
noOfTriangles += 4;
td1 = getTET3();
td2 = getTET3();
td3 = getTET3();
td4 = getTET3();
noOfTetrahedra += 3;
// Test which edge is refined
refEdges = 0;
ed = e1;
if ( (ed->firstSon)!=0 )
{
refEdges++;
p1=td->p1; p2=td->p2; p3=td->p3; p4=td->p4;
pm1 = e1->pm; pm2 = e5->pm;
t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
n1 = td->n1; n2 = td->n2; n3 = td->n3; n4 = td->n4;
}
ed = e2;
if ( (ed->firstSon)!=0 )
{
refEdges++;
if ( refEdges > 1) cout << "Alarm in MkGreen2AT3 \n";
p1=td->p3; p2=td->p1; p3=td->p2; p4=td->p4;
pm1 = e2->pm; pm2 = e6->pm;
t1 = trs[3]; t2 = trs[1]; t3 = trs[2]; t4 = trs[4];
n1 = td->n3; n2 = td->n1; n3 = td->n2; n4 = td->n4;
}
ed = e3;
if ( (ed->firstSon)!=0 )
{
refEdges++;
if ( refEdges > 1) cout << "Alarm in MkGreen2AT3 \n";
p1=td->p2; p2=td->p3; p3=td->p1; p4=td->p4;
pm1 = e3->pm; pm2 = e4->pm;
t1 = trs[2]; t2 = trs[3]; t3 = trs[1]; t4 = trs[4];
n1 = td->n2; n2 = td->n3; n3 = td->n1; n4 = td->n4;
}
// Set the fields for the new objects
MkGreen(t1);
t = t1->firstSon;
if ( PinTr(p4,t) )
{ t11 = t; t12 = t->next; }
else
{ t12 = t; t11 = t->next; }
MkGreen(t2);
t = t2->firstSon;
if ( PinTr(p1,t) )
{ t21 = t; t22 = t->next; }
else
{ t22 = t; t21 = t->next; }
MkGreen(t3);
t = t3->firstSon;
if ( PinTr(p2,t) )
{ t31 = t; t32 = t->next; }
else
{ t32 = t; t31 = t->next; }
MkGreen(t4);
t = t4->firstSon;
if ( PinTr(p3,t) )
{ t41 = t; t42 = t->next; }
else
{ t42 = t; t41 = t->next; }
edNew->p1 = pm1; edNew->p2 = pm2;
edNew->pm = 0;
edNew->type = T_GREEN_2A;
e1New = EdinTr(t12,p2,pm2);
e2New = EdinTr(t22,p4,pm1);
e3New = EdinTr(t11,p4,pm2);
e4New = EdinTr(t41,p2,pm1);
e5New = EdinTr(t11,p3,pm2);
e6New = EdinTr(t22,p3,pm1);
e7New = EdinTr(t32,p1,pm2);
e8New = EdinTr(t21,p1,pm1);
t1New->e1 = edNew; t1New->e2 = e2New; t1New->e3 = e3New;
SetTP(t1New,0); t1New->type = T_GREEN_2A;
t2New->e1 = e1New; t2New->e2 = e4New; t2New->e3 = edNew;
SetTP(t2New,0); t2New->type = T_GREEN_2A;
t3New->e1 = edNew; t3New->e2 = e5New; t3New->e3 = e6New;
SetTP(t3New,0); t3New->type = T_GREEN_2A;
t4New->e1 = edNew; t4New->e2 = e7New; t4New->e3 = e8New;
SetTP(t4New,0); t4New->type = T_GREEN_2A;
SetTdP(td1,t11,t22,t1New,t3New,td);
SetTdE(td1,t11,t22,t1New,t3New);
neighbourAtFace(td1,n1,t11);
neighbourAtFace(td1,n2,t22);
td1->type = T_GREEN_2A;
SetTdP(td2,t12,t2New,t41,t3New,td);
SetTdE(td2,t12,t2New,t41,t3New);
neighbourAtFace(td2,n1,t12);
neighbourAtFace(td2,n4,t41);
td2->type = T_GREEN_2A;
SetTdP(td3,t32,t21,t1New,t4New,td);
SetTdE(td3,t32,t21,t1New,t4New);
neighbourAtFace(td3,n3,t32);
neighbourAtFace(td3,n2,t21);
td3->type = T_GREEN_2A;
SetTdP(td4,t31,t42,t2New,t4New,td);
SetTdE(td4,t31,t42,t2New,t4New);
neighbourAtFace(td4,n4,t42);
neighbourAtFace(td4,n3,t31);
td4->type = T_GREEN_2A;
innerNeighbourAtFace(td1,t1New,td3);
innerNeighbourAtFace(td1,t3New,td2);
innerNeighbourAtFace(td2,t2New,td4);
innerNeighbourAtFace(td2,t3New,td1);
innerNeighbourAtFace(td3,t1New,td1);
innerNeighbourAtFace(td3,t4New,td4);
innerNeighbourAtFace(td4,t2New,td2);
innerNeighbourAtFace(td4,t4New,td3);
td->firstSon = td1;
td->type = T_GREEN_2A;
tetList[tetList.h]->add(td1);
tetList[tetList.h]->add(td2);
tetList[tetList.h]->add(td3);
tetList[tetList.h]->add(td4);
// destruction of the faces
delete t1New;
delete t2New;
delete t3New;
delete t4New;
if ( !(t1->onBoundary()) ) {
delete t11;
delete t12;
}
if ( !(t2->onBoundary()) ) {
delete t21;
delete t22;
}
if ( !(t3->onBoundary()) ) {
delete t31;
delete t32;
}
if ( !(t4->onBoundary()) ) {
delete t41;
delete t42;
}
td->returnTriangles(trs);
return True;
}
//-------------------------------------------------------------------------
int MESH3:: MkGreen2BTd(TET3 *td)
{
PT3 *p1=0, *p2=0, *p3=0, *p4=0;
TR3 *t=0, *t1=0, *t2=0, *t3=0, *t4=0;
TR3 *trs[1+4];
TR3 *t1New=0, *t2New=0;
TR3 *t11=0, *t12=0, *t21=0, *t22=0, *t23=0, *t31=0, *t32=0;
TET3 *td1=0, *td2=0, *td3=0;
PATCH *n1=0, *n2=0, *n3=0, *n4=0;
EDG3 *e1=0, *e2=0, *e3=0, *e4=0, *e5=0, *e6=0;
EDG3 *e1New=0, *e2New=0, *e3New=0, *e4New=0;
int i, depth = (td->depth)+1;
// Refinement of td into three green tetrahedra
if ( (td->mark == T_GREEN_2B) ) td->mark = False;
else return True;
td->getTriangles(trs);
// Get two new triangles and three tetrahedra,
// if not enough space: return False
t1New = new TR3;
t2New = new TR3;
noOfTriangles += 2;
t1New->depth = t2New->depth = depth;
td1 = getTET3();
td2 = getTET3();
td3 = getTET3();
noOfTetrahedra += 2;
// Testing which face is green-2B refined
// The actual case is transformed on the td with e4, e6 refined.
for (i=1; i<=4; i++)
{
t = trs[i];
if ( TestFace(t) > 1 )
{
MkGreenIrr(t);
if ( i==1 )
{
if ( (td->e2)->firstSon == 0)
{
t1 = trs[3]; t2 = trs[1]; t3 = trs[2]; t4 = trs[4];
n1 = td->n3; n2 = td->n1; n3 = td->n2; n4 = td->n4;
e1 = td->e2; e2 = td->e3; e3 = td->e1;
e4 = td->e5; e5 = td->e6; e6 = td->e4;
p1 = td->p3; p2 = td->p1; p3 = td->p2; p4 = td->p4;
}
else
if ( (td->e4)->firstSon == 0)
{
t1 = trs[4]; t2 = trs[1]; t3 = trs[3]; t4 = trs[2];
n1 = td->n4; n2 = td->n1; n3 = td->n3; n4 = td->n2;
e1 = td->e4; e2 = td->e1; e3 = td->e6;
e4 = td->e2; e5 = td->e3; e6 = td->e5;
p1 = td->p4; p2 = td->p1; p3 = td->p3; p4 = td->p2;
}
else
if ( (td->e5)->firstSon == 0)
{
t1 = trs[2]; t2 = trs[1]; t3 = trs[4]; t4 = trs[3];
n1 = td->n2; n2 = td->n1; n3 = td->n4; n4 = td->n3;
e1 = td->e5; e2 = td->e6; e3 = td->e3;
e4 = td->e4; e5 = td->e1; e6 = td->e2;
p1 = td->p2; p2 = td->p1; p3 = td->p4; p4 = td->p3;
}
}
else
if ( i==2 )
{
if ( (td->e1)->firstSon == 0)
{
t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
n1 = td->n1; n2 = td->n2; n3 = td->n3; n4 = td->n4;
e1 = td->e1; e2 = td->e2; e3 = td->e3;
e4 = td->e4; e5 = td->e5; e6 = td->e6;
p1 = td->p1; p2 = td->p2; p3 = td->p3; p4 = td->p4;
}
else
if ( (td->e4)->firstSon == 0)
{
t1 = trs[3]; t2 = trs[2]; t3 = trs[4]; t4 = trs[1];
n1 = td->n3; n2 = td->n2; n3 = td->n4; n4 = td->n1;
e1 = td->e4; e2 = td->e5; e3 = td->e2;
e4 = td->e6; e5 = td->e3; e6 = td->e1;
p1 = td->p3; p2 = td->p2; p3 = td->p4; p4 = td->p1;
}
else
if ( (td->e6)->firstSon == 0)
{
t1 = trs[4]; t2 = trs[2]; t3 = trs[1]; t4 = trs[3];
n1 = td->n4; n2 = td->n2; n3 = td->n1; n4 = td->n3;
e1 = td->e6; e2 = td->e3; e3 = td->e5;
e4 = td->e1; e5 = td->e2; e6 = td->e4;
p1 = td->p4; p2 = td->p2; p3 = td->p1; p4 = td->p3;
}
}
else
if ( i==3 )
{
if ( (td->e3)->firstSon == 0)
{
t1 = trs[2]; t2 = trs[3]; t3 = trs[1]; t4 = trs[4];
n1 = td->n2; n2 = td->n3; n3 = td->n1; n4 = td->n4;
e1 = td->e3; e2 = td->e1; e3 = td->e2;
e4 = td->e6; e5 = td->e4; e6 = td->e5;
p1 = td->p2; p2 = td->p3; p3 = td->p1; p4 = td->p4;
}
else
if ( (td->e5)->firstSon == 0)
{
t1 = trs[4]; t2 = trs[3]; t3 = trs[2]; t4 = trs[1];
n1 = td->n4; n2 = td->n3; n3 = td->n2; n4 = td->n1;
e1 = td->e5; e2 = td->e2; e3 = td->e4;
e4 = td->e3; e5 = td->e1; e6 = td->e6;
p1 = td->p4; p2 = td->p3; p3 = td->p2; p4 = td->p1;
}
else
if ( (td->e6)->firstSon == 0)
{
t1 = trs[1]; t2 = trs[3]; t3 = trs[4]; t4 = trs[2];
n1 = td->n1; n2 = td->n3; n3 = td->n4; n4 = td->n2;
e1 = td->e6; e2 = td->e4; e3 = td->e1;
e4 = td->e5; e5 = td->e2; e6 = td->e3;
p1 = td->p1; p2 = td->p3; p3 = td->p4; p4 = td->p2;
}
}
else
if ( i==4 )
{
if ( (td->e1)->firstSon == 0)
{
t1 = trs[3]; t2 = trs[4]; t3 = trs[1]; t4 = trs[2];
n1 = td->n3; n2 = td->n4; n3 = td->n1; n4 = td->n2;
e1 = td->e1; e2 = td->e6; e3 = td->e4;
e4 = td->e3; e5 = td->e5; e6 = td->e2;
p1 = td->p3; p2 = td->p4; p3 = td->p1; p4 = td->p2;
}
else
if ( (td->e2)->firstSon == 0)
{
t1 = trs[2]; t2 = trs[4]; t3 = trs[3]; t4 = trs[1];
n1 = td->n2; n2 = td->n4; n3 = td->n3; n4 = td->n1;
e1 = td->e2; e2 = td->e4; e3 = td->e5;
e4 = td->e1; e5 = td->e6; e6 = td->e3;
p1 = td->p2; p2 = td->p4; p3 = td->p3; p4 = td->p1;
}
else
if ( (td->e3)->firstSon == 0)
{
t1 = trs[1]; t2 = trs[4]; t3 = trs[2]; t4 = trs[3];
n1 = td->n1; n2 = td->n4; n3 = td->n2; n4 = td->n3;
e1 = td->e3; e2 = td->e5; e3 = td->e6;
e4 = td->e2; e5 = td->e4; e6 = td->e1;
p1 = td->p1; p2 = td->p4; p3 = td->p2; p4 = td->p3;
}
}
break;
}
}
// After Transform we can work on the reference td with refined e4, e6
MkGreen(t1);
t = t1->firstSon;
if ( PinTr(p4,t) )
{ t11 = t; t12 = t->next; }
else
{ t12 = t; t11 = t->next; }
MkGreen(t3);
t = t3->firstSon;
if ( PinTr(p4,t) )
{ t31 = t; t32 = t->next; }
else
{ t32 = t; t31 = t->next; }
// subfaces of the grenn-IIB refined face
t21 = t2->firstSon;
t23 = t21->next;
t22 = t23->next;
// Set the fields for the new objects
e1New = EdinTr(t12,p2,e4->pm);
e3New = EdinTr(t32,p2,e6->pm);
e2New = EdinTr(t21,e4->pm,e6->pm);
t1New->e1 = e1New; t1New->e2 = e2New; t1New->e3 = e3New;
SetTP(t1New,0); t1New->type = T_GREEN_2B;
SetTdP(td1,t11,t21,t1New,t31,td);
SetTdE(td1,t11,t21,t1New,t31);
neighbourAtFace(td1,n1,t11);
neighbourAtFace(td1,n2,t21);
neighbourAtFace(td1,n3,t31);
td1->type = T_GREEN_2B;
if ( PinTr(e6->pm,t23) )
{
e4New = EdinTr(t23,p3,e6->pm);
t2New->e1 = e2; t2New->e2 = e4New; t2New->e3 = e3New;
SetTP(t2New,0); t2New->type = T_GREEN_2B;
SetTdP(td2,t1New,t2New,t12,t22,td);
SetTdE(td2,t1New,t2New,t12,t22);
neighbourAtFace(td2,n1,t12);
neighbourAtFace(td2,n2,t22);
td2->type = T_GREEN_2B;
SetTdP(td3,t23,t32,t2New,t4,td);
SetTdE(td3,t23,t32,t2New,t4);
neighbourAtFace(td3,n2,t23);
neighbourAtFace(td3,n3,t32);
neighbourAtFace(td3,n4,t4);
td3->type = T_GREEN_2B;
}
else
{
e4New = EdinTr(t23,p1,e4->pm);
t2New->e1 = e1New; t2New->e2 = e4New; t2New->e3 = e3;
SetTP(t2New,0); t2New->type = T_GREEN_2B;
SetTdP(td2,t1New,t2New,t32,t22,td);
SetTdE(td2,t1New,t2New,t32,t22);
neighbourAtFace(td2,n3,t32);
neighbourAtFace(td2,n2,t22);
td2->type = T_GREEN_2B;
SetTdP(td3,t23,t12,t2New,t4,td);
SetTdE(td3,t23,t12,t2New,t4);
neighbourAtFace(td3,n2,t23);
neighbourAtFace(td3,n1,t12);
neighbourAtFace(td3,n4,t4);
td3->type = T_GREEN_2B;
}
innerNeighbourAtFace(td1,t1New,td2);
innerNeighbourAtFace(td2,t1New,td1);
innerNeighbourAtFace(td2,t2New,td3);
innerNeighbourAtFace(td3,t2New,td2);
td->firstSon = td1;
td->type = T_GREEN_2B;
tetList[tetList.h]->add(td1);
tetList[tetList.h]->add(td2);
tetList[tetList.h]->add(td3);
// destruction of the faces
delete t1New;
delete t2New;
if ( !(t1->onBoundary()) ) {
delete t11;
delete t12;
}
if ( !(t3->onBoundary()) ) {
delete t31;
delete t32;
}
if ( !(t2->onBoundary()) ) {
delete t21;
delete t22;
delete t23;
}
td->returnTriangles(trs);
return True;
}
//-------------------------------------------------------------------------
int MESH3:: MkGreen3Td(TET3 *td)
{
PT3 *p1, *p2, *p3, *p4;
PT3 *p6,*p7,*p9;
TR3 *t, *t1, *t2, *t3, *t4, *trs[1+4];
TR3 *s=0, *s1=0, *s2=0, *s3=0, *s4=0;
TR3 *t1New=0, *t2New=0, *t3New=0, *t4New=0, *t5New=0, *t6New=0, *t7New=0,
*t8New=0, *t9New=0;
TET3 *td1,*td2,*td3,*td4;
PATCH *n1, *n2, *n3, *n4;
EDG3 *e1=td->e1,*e2=td->e2,*e3=td->e3,*e4=td->e4, *e5=td->e5, *e6=td->e6;
EDG3 *e1New, *e2New,*e3New;
EDG3 *e14, *e24, *e34;
int i, depth = (td->depth)+1;
// Refinement of td into four green tetrahedra
if (td->mark == T_GREEN_3) td->mark = False;
else return True;
td->getTriangles(trs);
for (i=1; i<=4; i++)
{
t = trs[i];
if ( TestFace(t) > 1 )
{
RefineTr(t);
break;
}
}
// Get three new edges, nine new triangles and four tetrahedra,
// if not enough space: return False
t1New = new TR3;
t2New = new TR3;
t3New = new TR3;
noOfTriangles += 3;
t1New->depth = depth; t2New->depth = depth; t3New->depth = depth;
td1 = getTET3();
td2 = getTET3();
td3 = getTET3();
td4 = getTET3();
noOfTetrahedra += 3;
// Testing if one neighbor tetrahedron is refined and setting pi, ei
// corresponding to this neighbor
t = trs[1];
if ( (t->firstSon)!=0 && (t->type == T_RED) )
{
p1=td->p1; p2=td->p2; p3=td->p3; p4=td->p4;
p6 = e5->pm; p7 = e2->pm; p9 = e4->pm;
t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
n1 = td->n1; n2 = td->n2; n3 = td->n3; n4 = td->n4;
}
else
{
t = trs[2];
if ( (t->firstSon)!=0 && (t->type == T_RED) )
{
p1=td->p2; p2=td->p4; p3=td->p3; p4=td->p1;
p6 = e6->pm; p7 = e4->pm; p9 = e1->pm;
t1 = trs[2]; t2 = trs[4]; t3 = trs[3]; t4 = trs[1];
n1 = td->n2; n2 = td->n4; n3 = td->n3; n4 = td->n1;
}
else
{
t = trs[3];
if ( (t->firstSon)!=0 && (t->type == T_RED) )
{
p1=td->p3; p2=td->p4; p3=td->p1; p4=td->p2;
p6 = e5->pm; p7 = e6->pm; p9 = e3->pm;
t1 = trs[3]; t2 = trs[4]; t3 = trs[1]; t4 = trs[2];
n1 = td->n3; n2 = td->n4; n3 = td->n1; n4 = td->n2;
}
else
{
t = trs[4];
if ( (t->firstSon)!=0 && (t->type == T_RED) )
{
p1=td->p4; p2=td->p1; p3=td->p3; p4=td->p2;
p6 = e3->pm; p7 = e1->pm; p9 = e2->pm;
t1 = trs[4]; t2 = trs[1]; t3 = trs[3]; t4 = trs[2];
n1 = td->n4; n2 = td->n1; n3 = td->n3; n4 = td->n2;
}
else return False;
}
}
}
// Set the fields for the new objects
MkGreen(t2);
MkGreen(t3);
MkGreen(t4);
s = t1->firstSon;
for (i=0;i<4;i++)
{
if (PinTr(p2,s)) s2 = s;
else if (PinTr(p3,s)) s1 = s;
else if (PinTr(p4,s)) s3 = s;
else {s4 = s;}
s = s->next;
}
s = t2->firstSon;
e1New = EdinTr(s,p1,p9);
if (PinTr(p4,s))
{ t4New = s;
t5New = s->next; }
else
{ t5New = s;
t4New = s->next; }
s = t4->firstSon;
e2New = EdinTr(s,p1,p7);
if (PinTr(p2,s))
{ t7New = s;
t6New = s->next; }
else
{ t6New = s;
t7New = s->next; }
s = t3->firstSon;
e3New = EdinTr(s,p1,p6);
if (PinTr(p2,s))
{ t8New = s;
t9New = s->next; }
else
{ t9New = s;
t8New = s->next; }
e14 = EdinTr(s1,p7,p9);
e24 = EdinTr(s2,p6,p7);
e34 = EdinTr(s3,p6,p9);
t1New->e1 = e1New; t1New->e2 = e14; t1New->e3 = e2New;
t2New->e1 = e2New; t2New->e2 = e24; t2New->e3 = e3New;
t3New->e1 = e3New; t3New->e2 = e34; t3New->e3 = e1New;
SetTP(t1New,0); t1New->type = T_GREEN_3;
SetTP(t2New,0); t2New->type = T_GREEN_3;
SetTP(t3New,0); t3New->type = T_GREEN_3;
SetTdP(td1,s1,t1New,t5New,t6New,td);
SetTdE(td1,s1,t1New,t5New,t6New);
neighbourAtFace(td1,n1,s1);
neighbourAtFace(td1,n2,t5New);
neighbourAtFace(td1,n4,t6New);
td1->type = T_GREEN_3;
SetTdP(td2,s2,t2New,t7New,t8New,td);
SetTdE(td2,s2,t2New,t7New,t8New);
neighbourAtFace(td2,n1,s2);
neighbourAtFace(td2,n4,t7New);
neighbourAtFace(td2,n3,t8New);
td2->type = T_GREEN_3;
SetTdP(td3,s3,t3New,t9New,t4New,td);
SetTdE(td3,s3,t3New,t9New,t4New);
neighbourAtFace(td3,n1,s3);
neighbourAtFace(td3,n3,t9New);
neighbourAtFace(td3,n2,t4New);
td3->type = T_GREEN_3;
SetTdP(td4,s4,t1New,t2New,t3New,td);
SetTdE(td4,s4,t1New,t2New,t3New);
neighbourAtFace(td4,n1,s4);
td4->type = T_GREEN_3;
innerNeighbourAtFace(td1,t1New,td4);
innerNeighbourAtFace(td2,t2New,td4);
innerNeighbourAtFace(td3,t3New,td4);
innerNeighbourAtFace(td4,t1New,td1);
innerNeighbourAtFace(td4,t2New,td2);
innerNeighbourAtFace(td4,t3New,td3);
td->firstSon = td1;
td->type = T_GREEN_3;
tetList[tetList.h]->add(td1);
tetList[tetList.h]->add(td2);
tetList[tetList.h]->add(td3);
tetList[tetList.h]->add(td4);
// destruction of the faces
delete t1New;
delete t2New;
delete t3New;
if ( !(t1->onBoundary()) ) {
delete s1;
delete s2;
delete s3;
delete s4;
}
if ( !(t2->onBoundary()) ) {
delete t4New;
delete t5New;
}
if ( !(t3->onBoundary()) ) {
delete t8New;
delete t9New;
}
if ( !(t4->onBoundary()) ) {
delete t6New;
delete t7New;
}
td->returnTriangles(trs);
return True;
}
//-------------------------------------------------------------------------
void MESH3:: RelGreen()
{
TR3 *t1old, *t2old, *tFather;
EDG3 *eMid = 0, *ed, *edNext;
int k, depth;
// remove all green triangles (on the boundary)
FORALL_DOWNWARD(trList,depth)
{
TR3 *t = trList[depth]->first;
while (t!=0)
{
k = 0;
tFather = t->father;
if (tFather==0) { t = t->next; continue; }
if ((tFather->type)!=T_GREEN_1) { t = t->next; continue; }
// Remove two green triangles, print internal error message
// if the green edge is not found
t1old = t; t2old = t->next;
if (((t->e1)->type)==T_GREEN_1) { eMid = t->e1; k++;}
if (((t->e2)->type)==T_GREEN_1) { eMid = t->e2; k++;}
if (((t->e3)->type)==T_GREEN_1) { eMid = t->e3; k++;}
if (k!=1)
{
cout << "\n*** RemoveGreen: Triangle with impossible state k="
<< k << "\n";
}
t = t2old->next;
// Remove the two green triangles; return the space
trList[t1old->depth]->remove(t1old);
trList[t2old->depth]->remove(t2old);
delete t1old;
delete t2old;
noOfTriangles--;
tFather->firstSon = 0;
tFather->type = T_WHITE;
tFather->mark = False;
edgList[eMid->depth]->remove(eMid);
edgAlloc.Return(eMid);
noOfEdges--;
}
}
FORALL_DOWNWARD(edgList,depth)
{
ed = edgList[depth]->first;
while (ed!=0)
{
edNext = ed->next;
if ( (ed->type)==T_GREEN_1 || ((ed->type)==T_GREEN_3) )
{
edgList[ed->depth]->remove(ed);
edgAlloc.Return(ed);
noOfEdges--;
noOfTriangles--;
}
ed = edNext;
}
}
return;
}
//-------------------------------------------------------------------------
void MESH3:: RelGreenIrr()
{
TR3 *t1,*t2,*t3,*tFather;
EDG3 *e1Mid = 0, *e2Mid = 0, *ed, *edNext;
int k, depth;
// remove all irregular green triangles ( on the boundary)
FORALL_DOWNWARD(trList,depth)
{
TR3 *t = trList[depth]->first;
// -- Loop for all triangles; ApplyT could not be used
while (t != 0)
{
tFather = t->father;
if (tFather==0) { t = t->next; continue; }
if (tFather->type != T_GREEN_2B) { t = t->next; continue; }
// -- Remove three green triangles, print internal error message
// if the green edge is not found
t1 = t; t2 = t1->next; t3 = t2->next;
k = 0;
if (((t1->e1)->type)==T_GREEN_2B) { e1Mid = t1->e1; k++;}
if (((t1->e2)->type)==T_GREEN_2B) { e1Mid = t1->e2; k++;}
if (((t1->e3)->type)==T_GREEN_2B) { e1Mid = t1->e3; k++;}
if (k != 1)
cout << "\n*** 1RelGreenIrr: Triangle with impossible state k="
<< k << "\n";
k = 0;
if (((t2->e1)->type)==T_GREEN_2B) { e2Mid = t2->e1; k++;}
if (((t2->e2)->type)==T_GREEN_2B) { e2Mid = t2->e2; k++;}
if (((t2->e3)->type)==T_GREEN_2B) { e2Mid = t2->e3; k++;}
if (k!=1)
cout << "\n*** 2RelGreenIrr: Triangle with impossible state k="
<< k << "\n";
t = t3->next; // next while-step
// Remove the three green triangles; return the space
trList[t1->depth]->remove(t1);
trList[t2->depth]->remove(t2);
trList[t3->depth]->remove(t3);
delete t1;
delete t2;
delete t3;
noOfTriangles -= 2;
tFather->firstSon = 0;
tFather->type = T_WHITE;
edgList[e1Mid->depth]->remove(e1Mid);
edgList[e2Mid->depth]->remove(e2Mid);
edgAlloc.Return(e1Mid);
edgAlloc.Return(e2Mid);
noOfEdges -= 2;
}
}
FORALL_DOWNWARD(edgList,depth)
{
ed = edgList[depth]->first;
while (ed!=0)
{
edNext = ed->next;
if ((ed->type)!=T_GREEN_2B) { ed = edNext; continue; }
edgList[ed->depth]->remove(ed);
edgAlloc.Return(ed);
noOfEdges--;
noOfTriangles -= 1;
ed = edNext;
}
}
return;
}
//-------------------------------------------------------------------------
void MESH3:: RelGreen1Td()
{
TET3 *td1old,*td2old,*tdFather;
TET3 *td = tetList[refStep]->first;
// remove all Green-I tetrahedra
while (td != 0)
{
tdFather = td->father;
if (tdFather == 0) { td = td->next; continue; }
if (tdFather->type != T_GREEN_1) { td = td->next; continue; }
// -- two Green-I tetrahedra
td1old = td; td2old = td->next;
td = td2old->next;
if (td1old->mark) tdFather->mark = True; // refine father
else if (td2old->mark) tdFather->mark = True; // refine father
else tdFather->mark = False;
// -- correct the white neighbors
correctWhiteNbs(td1old);
correctWhiteNbs(td2old);
// -- Remove the two Green-I tetrahedra from the list of tetrahedra
tetList[refStep]->remove(td1old);
tetList[refStep]->remove(td2old);
returnTET3(td1old);
returnTET3(td2old);
tdFather->firstSon = 0;
tdFather->type = T_WHITE;
noOfTetrahedra--;
noOfTriangles--;
}
return;
}
//-------------------------------------------------------------------------
void MESH3:: RelGreen2ATd()
{
TET3 *td1old,*td2old,*td3old,*td4old,*tdFather;
EDG3 *edMid=0;
TET3 *td = tetList[refStep]->first;
// remove all green-2A tetrahedra
while (td != 0)
{
tdFather = td->father;
if (tdFather==0) { td = td->next; continue; }
if (tdFather->type != T_GREEN_2A) { td = td->next; continue; }
// Remove four green-2A tetrahedra, print internal error message
// if the inner green-2A edge is not found
td1old = td; td2old = td->next;
td3old = td2old->next; td4old = td3old->next;
td = td4old->next;
if (td1old->mark) tdFather->mark = True; // refine father
else if (td2old->mark) tdFather->mark = True; // refine father
else if (td3old->mark) tdFather->mark = True; // refine father
else if (td4old->mark) tdFather->mark = True; // refine father
else tdFather->mark = False;
// inner edge
if ((td1old->e1)->type == T_GREEN_2A) edMid = td1old->e1;
else
if ((td1old->e2)->type == T_GREEN_2A) edMid = td1old->e2;
else
if ((td1old->e3)->type == T_GREEN_2A) edMid = td1old->e3;
else
if ((td1old->e4)->type == T_GREEN_2A) edMid = td1old->e4;
else
if ((td1old->e5)->type == T_GREEN_2A) edMid = td1old->e5;
else
if ((td1old->e6)->type == T_GREEN_2A) edMid = td1old->e6;
if (edMid == 0) cout << "\n\nAlarm in RelGreen2ATd: edMid == 0\n";
// -- correct the white neighbours
correctWhiteNbs(td1old);
correctWhiteNbs(td2old);
correctWhiteNbs(td3old);
correctWhiteNbs(td4old);
// Remove the four green ones from the list of tetrahedra
tetList[refStep]->remove(td1old);
tetList[refStep]->remove(td2old);
tetList[refStep]->remove(td3old);
tetList[refStep]->remove(td4old);
returnTET3(td1old);
returnTET3(td2old);
returnTET3(td3old);
returnTET3(td4old);
tdFather->firstSon = 0;
tdFather->type = T_WHITE;
noOfTetrahedra -= 3;
noOfTriangles -= 4;
edgList[edMid->depth]->remove(edMid);
edgAlloc.Return(edMid);
noOfEdges -=1;
}
return;
}
//-------------------------------------------------------------------------
void MESH3:: RelGreen2BTd()
{
TET3 *td1old,*td2old,*td3old,*tdFather;
TET3 *td = tetList[refStep]->first;
// remove all green-2B tetrahedra
while (td != 0)
{
tdFather = td->father;
if (tdFather == 0) { td = td->next; continue; }
if (tdFather->type != T_GREEN_2B) { td = td->next; continue; }
// Remove three green2B tetrahedra
td1old = td; td2old = td->next; td3old = td2old->next;
td = td3old->next; // next while-step
if (td1old->mark) tdFather->mark = True; // refine father
else if (td2old->mark) tdFather->mark = True; // refine father
else if (td3old->mark) tdFather->mark = True; // refine father
else tdFather->mark = False;
// -- correct the white neighbours
correctWhiteNbs(td1old);
correctWhiteNbs(td2old);
correctWhiteNbs(td3old);
// Remove the three green ones from the list of tetrahedra
tetList[refStep]->remove(td1old);
tetList[refStep]->remove(td2old);
tetList[refStep]->remove(td3old);
returnTET3(td1old);
returnTET3(td2old);
returnTET3(td3old);
tdFather->firstSon = 0;
tdFather->type = T_WHITE;
noOfTetrahedra -=2;
noOfTriangles -= 2;
}
return;
}
//-------------------------------------------------------------------------
void MESH3:: RelGreen3Td()
{
TET3 *td1old,*td2old,*td3old,*td4old,*tdFather;
TET3 *td = tetList[refStep]->first;
// remove all green-III tetrahedra
while (td != 0)
{
tdFather = td->father;
if (tdFather==0) { td = td->next; continue; }
if ((tdFather->type)!=T_GREEN_3) { td = td->next; continue; }
// Remove four green tetrahedra, print internal error message
// if the inner green faces are not found
td1old = td; td2old = td->next;
td3old = td2old->next; td4old = td3old->next;
td = td4old->next; // next step in while
if (td1old->mark) tdFather->mark = True; // refine father
else if (td2old->mark) tdFather->mark = True; // refine father
else if (td3old->mark) tdFather->mark = True; // refine father
else if (td4old->mark) tdFather->mark = True; // refine father
else tdFather->mark = False;
// -- correct the white neighbours
correctWhiteNbs(td1old);
correctWhiteNbs(td2old);
correctWhiteNbs(td3old);
correctWhiteNbs(td4old);
// Substitute the old (father) tetrahedron for the four green ones
// in the list of tetrahedra; return the space
tetList[refStep]->remove(td1old);
tetList[refStep]->remove(td2old);
tetList[refStep]->remove(td3old);
tetList[refStep]->remove(td4old);
returnTET3(td1old);
returnTET3(td2old);
returnTET3(td3old);
returnTET3(td4old);
tdFather->firstSon = 0;
tdFather->type = T_WHITE;
noOfTetrahedra -= 3;
noOfTriangles -= 3;
}
return;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool TET3:: noTdinNb(TET3* tet)
{
TET* nbTet;
nbTet = n1->castToTET3();
if ( nbTet == tet) return True;
nbTet = n2->castToTET3();
if ( nbTet == tet) return True;
nbTet = n3->castToTET3();
if ( nbTet == tet) return True;
nbTet = n4->castToTET3();
if ( nbTet == tet) return True;
return False;
}
//-------------------------------------------------------------------------
void TET3:: correctNbinTet(TET3 *tet)
{
TET* nbTet;
nbTet = n1->castToTET3();
if ( nbTet == tet) {n1 = tet->father; return; }
nbTet = n2->castToTET3();
if ( nbTet == tet) {n2 = tet->father; return; }
nbTet = n3->castToTET3();
if ( nbTet == tet) {n3 = tet->father; return; }
nbTet = n4->castToTET3();
if ( nbTet == tet) {n4 = tet->father; return; }
// cout << "\nImpossible state in correctNbinTet\n";
}
void MESH3:: correctWhiteNbs(TET3 *tet)
{
TET3* nbTet;
Vector<PATCH*> nbs(4);
int i;
tet->getNeighbours(nbs);
FORALL(nbs,i)
{
if (nbs[i])
{
nbTet = nbs[i]->castToTET3();
nbTet->correctNbinTet(tet);
// if ( (nbTet->type != T_GREEN_1) && (nbTet->type != T_GREEN_2A) &&
// (nbTet->type != T_GREEN_2B) && (nbTet->type != T_GREEN_3) )
}
}
}
//-------------------------------------------------------------------------
int MESH3:: CheckTd(TET3 *td)
{
TET3 *nbTd;
PATCH *n1=td->n1, *n2=td->n2, *n3=td->n3, *n4=td->n4;
if ( n1 == 0 )
{
cout << "\n\tneighbour 1==0:\n";
td->print(cout);
(td->p2)->print(cout); (td->p3)->print(cout); (td->p4)->print(cout);
}
if ( n2 == 0 )
{
cout << "\n\tneighbour 2==0:\n";
td->print(cout);
(td->p1)->print(cout); (td->p3)->print(cout); (td->p4)->print(cout);
}
if ( n3 == 0 )
{
cout << "\n\tneighbour 3==0:\n";
td->print(cout);
(td->p1)->print(cout); (td->p2)->print(cout); (td->p4)->print(cout);
}
if ( n4 == 0 )
{
cout << "\n\tneighbour 4==0:\n";
td->print(cout);
(td->p1)->print(cout); (td->p2)->print(cout); (td->p3)->print(cout);
}
nbTd = n1->castToTET3();
if ( nbTd )
{
if( !nbTd->noTdinNb(td) ) {
cout << "\n\t error in neighbor relation\n";
td->print(cout); nbTd->print(cout);
}
}
nbTd = n2->castToTET3();
if ( nbTd )
{
if( !nbTd->noTdinNb(td) ) {
cout << "\n\t error in neighbor relation\n";
td->print(cout); nbTd->print(cout);
}
}
nbTd = n3->castToTET3();
if ( nbTd )
{
if( !nbTd->noTdinNb(td) ) {
cout << "\n\t error in neighbor relation\n";
td->print(cout); nbTd->print(cout);
}
}
nbTd = n4->castToTET3();
if ( nbTd )
{
if( !nbTd->noTdinNb(td) ) {
cout << "\n\t error in neighbor relation\n";
td->print(cout); nbTd->print(cout);
}
}
if ( n1->castToTET3() && n2->castToTET3() &&
n3->castToTET3() && n4->castToTET3() )
{
if (td->boundP != INTERIOR) {
cout << "\n\t error in boundary definition: boundP != INTERIOR\n";
td->print(cout);
}
}
else
{
if (td->boundP != BOUNDARY) {
cout << "\n\t error in boundary definition: boundP != BOUNDARY\n";
td->print(cout);
}
}
return True;
}
//-------------------------------------------------------------------------
int MESH3:: consistencyCheck()
{
PT3 *p;
TET3 *tet;
//int no, noInner, noBound, noWithFather, noWithoutFather;
//int noWhite, noRed, noGreenI, noGreenIIA, noGreenIIB, noGreenIII;
DListIter<TET3> iterTET(tetras());
cout << "\n\t Begin of consistencyCheck:\n";
cout << "\n\t no of points = " << noOfPoints << "\n";
cout << "\t no of edges = " << noOfEdges << "\n";
cout << "\t no of triangles = " << noOfTriangles << "\n";
cout << "\t no of tetrahedra = " << noOfTetrahedra << "\n";
int eulerC = noOfPoints - noOfEdges + noOfTriangles - noOfTetrahedra;
cout << "\n\t Euler's characteristic = " << eulerC << "\n";
if (eulerC != 1) cout << "\n\t Dangerous: Euler's characteristic not 1\n\n";
Vector<int> PNodes(noOfPoints);
// saving the node numbers for later reuse
int i = 0;
DListIter<PT3> PIter(points());
while ((p=PIter.all()))
{
PNodes[++i] = p->node;
p->node = i;
}
while ((tet=iterTET.all())) {
if(!CheckTd(tet)) return False;;
}
// restore of node numbers
i = 0;
PIter.reset();
while ((p=PIter.all())) p->node = PNodes[++i];
CompAngles();
/*
no = 0; noInner = 0; noBound = 0;
while (p =iterPT.all()) {
no++;
if ( p->boundP == DIRICHLET ) noBound++;
if ( p->boundP == INTERIOR ) noInner++;
}
cout << "\n\t no of points in list = " << no;
cout << "\n\t no of boundary points = " << noBound;
cout << "\n\t no of inner points = " << noInner;
no = 0; noInner = 0; noBound = 0; noWithFather=0; noWithoutFather=0;
noWhite=0; noRed=0; noGreenI=0; noGreenIIA=0; noGreenIIB=0; noGreenIII=0;
noZero=0; noGreen=0;
while (e =iterEDG.all()) {
no++;
if ( e->boundP == DIRICHLET ) noBound++;
if ( e->boundP == INTERIOR ) noInner++;
if ( e->type == T_WHITE ) noWhite++;
if ( e->type == T_RED ) noRed++;
if ( e->type == T_GREEN_1 ) noGreenI++;
if ( e->type == T_GREEN_2A) noGreenIIA++;
if ( e->type == T_GREEN_2B) noGreenIIB++;
if ( e->type == T_GREEN_3 ) noGreenIII++;
if ( e->type == T_GREEN ) noGreen++;
if ( e->type == T_ZERO ) noZero++;
if ( e->father != 0 ) noWithFather++;
else if ( e->father == 0 ) noWithoutFather++;
}
cout << "\n\t no of edges in list = " << no;
cout << "\n\t no of boundary edges = " << noBound;
cout << "\n\t no of inner edges = " << noInner;
cout << "\n\t no of edges with father = " << noWithFather;
cout << "\n\t no of edges without father = " << noWithoutFather;
cout << "\n\t no of white edges = " << noWhite;
cout << "\n\t no of red edges = " << noRed;
cout << "\n\t no of greenI edges = " << noGreenI;
cout << "\n\t no of greenIIA edges = " << noGreenIIA;
cout << "\n\t no of greenIIB edges = " << noGreenIIB;
cout << "\n\t no of greenIII edges = " << noGreenIII;
cout << "\n\t no of green edges = " << noGreen;
cout << "\n\t no of zero edges = " << noZero;
no = 0;
iterTET.reset();
while (tet=iterTET.all()) no++;
cout << "\n\t no of tets in list = " << no << "\n";
*/
cout << "\n\t End of consistencyCheck.\n";
cout << "\n\t ----------------------------------------------------------\n\n";
return True;
}
//-------------------------------------------------------------------------
int MESH3:: CheckAngle(TR3 *t)
{
PT3 *p1=t->p1, *p2=t->p2, *p3=t->p3;
Real e1X, e1Y, e1Z;
Real e2X, e2Y, e2Z;
Real s12, s11, s22;
Real alpha1, alpha2, alpha3;
int classI;
e1X = p2->x - p1->x;
e1Y = p2->y - p1->y;
e1Z = p2->z - p1->z;
e2X = p3->x - p1->x;
e2Y = p3->y - p1->y;
e2Z = p3->z - p1->z;
s12 = e1X*e2X + e1Y*e2Y + e1Z*e2Z;
s11 = sqrt(e1X*e1X + e1Y*e1Y + e1Z*e1Z);
s22 = sqrt(e2X*e2X + e2Y*e2Y + e2Z*e2Z);
alpha1 = 180.0*acos(s12/(s11*s22))/REALPI;
if (alpha1 > alphaMax) alphaMax = alpha1;
if (alpha1 < alphaMin) alphaMin = alpha1;
e1X = p1->x - p2->x;
e1Y = p1->y - p2->y;
e1Z = p1->z - p2->z;
e2X = p3->x - p2->x;
e2Y = p3->y - p2->y;
e2Z = p3->z - p2->z;
s12 = e1X*e2X + e1Y*e2Y + e1Z*e2Z;
s11 = sqrt(e1X*e1X + e1Y*e1Y + e1Z*e1Z);
s22 = sqrt(e2X*e2X + e2Y*e2Y + e2Z*e2Z);
alpha2 = 180.0*acos(s12/(s11*s22))/REALPI;
if (alpha2 > alphaMax) alphaMax = alpha2;
if (alpha2 < alphaMin) alphaMin = alpha2;
alpha3 = 180.0-alpha1-alpha2;
if (alpha3 > alphaMax) alphaMax = alpha3;
if (alpha3 < alphaMin) alphaMin = alpha3;
classI = (int)(alpha1/5);
alpha[classI]++;
classI = (int)(alpha2/5);
alpha[classI]++;
classI = (int)(alpha3/5);
alpha[classI]++;
if ( (alpha1 < 10.) || (alpha2 < 10.) || (alpha3 < 10.) )
{
cout << " alphas: " << alpha1 << alpha2 << alpha3 << "\n";
cout << " P1 : " << p1->x << p1->y << p1->z << "\n";
cout << " P2 : " << p2->x << p2->y << p2->z << "\n";
cout << " P3 : " << p3->x << p3->y << p3->z << "\n";
}
return True;
}
//-------------------------------------------------------------------------
int MESH3:: CheckCircle(TET3 *tet)
{
PT3 *p1, *p2, *p3, *p4;
Real e1X, e1Y, e1Z;
Real e2X, e2Y, e2Z;
Real dyz, dxz, dxy;
Real len, rhoT, sigma, vT, sT, hT;
EDG3 *eds[1+6];
int i;
TR3 *t1, *t2, *t3, *t4;
TR3 *trs[1+4];
Real x1, x2, x3, y1, y2, y3, z1, z2, z3;
tet->getTriangles(trs);
t1 = trs[1]; t2 = trs[2]; t3 = trs[3]; t4 = trs[4];
p1 = tet->p1;
p2 = tet->p2;
p3 = tet->p3;
p4 = tet->p4;
x1 = (p2->x) - (p1->x);
x2 = (p2->y) - (p1->y);
x3 = (p2->z) - (p1->z);
y1 = (p3->x) - (p1->x);
y2 = (p3->y) - (p1->y);
y3 = (p3->z) - (p1->z);
z1 = (p4->x) - (p1->x);
z2 = (p4->y) - (p1->y);
z3 = (p4->z) - (p1->z);
vT = SIXTH*Abs(x1*(y2*z3 - y3*z2) - x2*(y1*z3 -y3*z1) + x3*(y1*z2 - y2*z1));
// = tet->volume;
sT = ZERO;
p1=t1->p1; p2=t1->p2, p3=t1->p3;
e1X = p2->x - p1->x;
e1Y = p2->y - p1->y;
e1Z = p2->z - p1->z;
e2X = p3->x - p1->x;
e2Y = p3->y - p1->y;
e2Z = p3->z - p1->z;
dyz = e1Y*e2Z - e1Z*e2Y;
dxz = e1Z*e2X - e1X*e2Z;
dxy = e1X*e2Y - e1Y*e2X;
sT += sqrt(dyz*dyz + dxz*dxz + dxy*dxy);
p1=t2->p1; p2=t2->p2, p3=t2->p3;
e1X = p2->x - p1->x;
e1Y = p2->y - p1->y;
e1Z = p2->z - p1->z;
e2X = p3->x - p1->x;
e2Y = p3->y - p1->y;
e2Z = p3->z - p1->z;
dyz = e1Y*e2Z - e1Z*e2Y;
dxz = e1Z*e2X - e1X*e2Z;
dxy = e1X*e2Y - e1Y*e2X;
sT += sqrt(dyz*dyz + dxz*dxz + dxy*dxy);
p1=t3->p1; p2=t3->p2, p3=t3->p3;
e1X = p2->x - p1->x;
e1Y = p2->y - p1->y;
e1Z = p2->z - p1->z;
e2X = p3->x - p1->x;
e2Y = p3->y - p1->y;
e2Z = p3->z - p1->z;
dyz = e1Y*e2Z - e1Z*e2Y;
dxz = e1Z*e2X - e1X*e2Z;
dxy = e1X*e2Y - e1Y*e2X;
sT += sqrt(dyz*dyz + dxz*dxz + dxy*dxy);
p1=t4->p1; p2=t4->p2, p3=t4->p3;
e1X = p2->x - p1->x;
e1Y = p2->y - p1->y;
e1Z = p2->z - p1->z;
e2X = p3->x - p1->x;
e2Y = p3->y - p1->y;
e2Z = p3->z - p1->z;
dyz = e1Y*e2Z - e1Z*e2Y;
dxz = e1Z*e2X - e1X*e2Z;
dxy = e1X*e2Y - e1Y*e2X;
sT += sqrt(dyz*dyz + dxz*dxz + dxy*dxy);
sT = 0.5*sT;
rhoT =6.0*vT/sT;
hT =ZERO;
tet->getEdges(eds);
for (i=1;i<=6;i++)
{
e1X = (eds[i]->p1)->x;
e1Y = (eds[i]->p1)->y;
e1Z = (eds[i]->p1)->z;
e2X = (eds[i]->p2)->x;
e2Y = (eds[i]->p2)->y;
e2Z = (eds[i]->p2)->z;
e2X = e2X - e1X;
e2Y = e2Y - e1Y;
e2Z = e2Z - e1Z;
len = sqrt(e2X*e2X + e2Y*e2Y + e2Z*e2Z);
if (len > hT) hT = len;
}
sigma = hT/rhoT;
if (sigma > sigmaMax) sigmaMax = sigma;
tet->returnTriangles(trs);
return True;
}
//-------------------------------------------------------------------------
int MESH3:: CompAngles()
{
TET3 *tet;
TR3 *t;
int i, noOfAngles=0;
for (i=0;i<36;i++)
{
alpha[i] = 0;
}
alphaMax = ZERO;
alphaMin = 90.0;
DListIter<TR3> iterTR(triangles());
while ((t=iterTR.all())) CheckAngle(t);
cout << "\n\t ----------------------------------------------------------\n";
cout << "\t angles on boundary faces: \n";
cout << "\n\t max. angle = " << alphaMax << "\n";
cout << "\t min. angle = " << alphaMin << "\n\n";
for (i=0;i<36;i++)
{
cout << "\t no of angles between " << i*5 << "and " << (i+1)*5 << ": "
<< alpha[i] << "\n";
noOfAngles += alpha[i];
}
cout << "\n\t total no of angles = " << noOfAngles << "\n";
cout << "\n\t ----------------------------------------------------------\n";
sigmaMax = ZERO;
DListIter<TET3> iterTET(tetras());
while ((tet=iterTET.all())) CheckCircle(tet);
cout << "\t sigmaT = " << sigmaMax << "\n";
cout << "\n\t ----------------------------------------------------------\n";
return True;
}
//-------------------------------------------------------------------------
Real MESH3:: Length(PT3 *p1, PT3 *p2)
{
double sum = ZERO, diff;
diff = (p1->x) - (p2->x);
sum += diff * diff;
diff = (p1->y) - (p2->y);
sum += diff * diff;
diff = (p1->z) - (p2->z);
sum += diff * diff;
return sum;
}
//-------------------------------------------------------------------------
int MESH3:: PinTr(PT3 *p, TR3 *tr)
{
PT3 *p1 = tr->p1, *p2 = tr->p2, *p3 = tr->p3;
if (p==p1 || p==p2 || p==p3) return True;
else return False;
}
//-------------------------------------------------------------------------
int MESH3:: PinTet(PT3 *p, TET3 *tet)
{
PT3 *p1 = tet->p1, *p2 = tet->p2, *p3 = tet->p3, *p4 = tet->p4;
if (p==p1 || p==p2 || p==p3 || p==p4) return True;
else return False;
}
//-------------------------------------------------------------------------
EDG3* MESH3:: EdinTr(TR3 *tr,PT3 *p1,PT3 *p2)
{
EDG3 *e1 = tr->e1, *e2 = tr->e2, *e3 = tr->e3;
PT3 *pt1, *pt2;
pt1 = e1->p1; pt2 = e1->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) ) return e1;
pt1 = e2->p1; pt2 = e2->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) ) return e2;
pt1 = e3->p1; pt2 = e3->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) ) return e3;
cout << "\nEdinTr(...): edge not found\n";
cout.flush(); abort();
return 0;
}
//-------------------------------------------------------------------------
EDG3* MESH3:: EdinTet(PT3 *p1,PT3 *p2,TET3 *td)
{
EDG3 *e1 = td->e1, *e2 = td->e2, *e3 = td->e3,
*e4 = td->e4, *e5 = td->e5, *e6 = td->e6;
PT3 *pt1, *pt2;
pt1 = e1->p1; pt2 = e1->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
return e1;
pt1 = e2->p1; pt2 = e2->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
return e2;
pt1 = e3->p1; pt2 = e3->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
return e3;
pt1 = e4->p1; pt2 = e4->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
return e4;
pt1 = e5->p1; pt2 = e5->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
return e5;
pt1 = e6->p1; pt2 = e6->p2;
if ( ((p1 == pt1)&&(p2 == pt2)) || ((p2 == pt1)&&(p1 == pt2)) )
return e6;
return 0;
}
//-------------------------------------------------------------------------
int MESH3:: TestFace(TR3 *tr)
{
EDG3 *eds[1+3];
int i, noOfRefEdg = 0;
tr->getEdges(eds);
for (i=1; i<=3; i++)
{
if ( (eds[i]->firstSon) != 0)
noOfRefEdg++;
}
return noOfRefEdg;
}
//-------------------------------------------------------------------------
int MESH3:: TestEdsInTd(TET3 *tet)
{
EDG3 *eds[1+6];
int i, noOfRefEdg = 0;
tet->getEdges(eds);
for (i=1; i<=6; i++)
{
if ( (eds[i]->firstSon) != 0)
noOfRefEdg++;
}
return noOfRefEdg;
}
//-------------------------------------------------------------------------
int MESH3:: testSons(TR3 *t, TET3 *tet)
{
EDG3 *e1=t->e1, *e2=t->e2, *e3=t->e3, *ed=0;
PT3 *pm1, *pm2, *pm3;
TET3 *tetSon;
int noOfRefEdg = 0, i;
pm1 = e1->pm; pm2 = e2->pm; pm3 = e3->pm;
tetSon = tet->firstSon;
for (i=1; i<=8; i++)
{
ed = EdinTet(pm1,pm2,tetSon);
if (ed != 0)
if (ed->refined()) noOfRefEdg++;
ed = EdinTet(pm1,pm3,tetSon);
if (ed != 0)
if (ed->refined()) noOfRefEdg++;
ed = EdinTet(pm2,pm3,tetSon);
if (ed != 0)
if (ed->refined()) noOfRefEdg++;
if ( noOfRefEdg > 0 ) return noOfRefEdg;
tetSon = tetSon->next;
}
return noOfRefEdg;
}
//-------------------------------------------------------------------------
ostream& operator << (ostream& os, MESH3& t)
{
PT3* p;
EDG3* e;
TR3* tr;
TET3 * tet;
os << "\nTriangulation 3D";
if (t.fileName) os << " " << t.fileName;
os << ":\n";
DListIter<PT3> iterPT (t.points());
DListIter<TR3> iterTR (t.triangles());
DListIter<EDG3> iterEDG(t.edges());
DListIter<TET3> iterTET(t.tetras());
while ((p =iterPT.all())) p->print(cout);
while ((e =iterEDG.all())) e->print(cout);
while ((tr=iterTR.all())) tr->print(cout);
while ((tet=iterTET.all())) tet->print(cout);
os << "\nRefinement step " << t.refStep;
os << " -- points: " << t.noOfPoints;
os << " edges: " << t.noOfEdges;
os << " triangles: " << t.noOfTriangles;
os << " tetras: " << t.noOfTetrahedra;
os << "\n";
return os;
}
//-------------------------------------------------------------------------
void MESH3:: Info()
{
cout << "\n Triangulation: points edges tri's tetra's "
<< " Depth=" << maxDepth << " Steps=" << refStep
<< Form("\n%15s%8d%8d%8d%8d", " ",
noOfPoints, noOfEdges, noOfTriangles, noOfTetrahedra)
<< "\n";
}
//-------------------------------------------------------------------------
int MESH3:: MemSpace()
{
int space = varAllocator.MemSpace();
space += tetAlloc.MemSpace();
space += edgAlloc.MemSpace();
space += ptAlloc.MemSpace();
space += trList[0]->first->alloc.MemSpace();
return space;
}
//-------------------------------------------------------------------------
int EFMESH3:: MemSpace()
{
int space = varAllocator.MemSpace();
space += EFTetAlloc.MemSpace();
space += edgAlloc.MemSpace();
space += ptAlloc.MemSpace();
space += trList[0]->first->alloc.MemSpace();
return space;
}
//-------------------------------------------------------------------------
void MESH3:: refinementInfo(int* noPt, int* noEdg, int* noTr, int* noTet)
{
int blocks = varAllocator.NoOfBlocks();
blocks += tetAlloc.NoOfBlocks();
blocks += edgAlloc.NoOfBlocks();
blocks += ptAlloc.NoOfBlocks();
blocks += trList[0]->first->alloc.NoOfBlocks();
int space = MemSpace();
cout << "\n Refine: points edges tri's tetra's "
<< " Depth=" << maxDepth << " Steps=" << refStep
<< Form("\n%15s%8d%8d%8d%8d", "previous",
noPt[0], noEdg[0], noTr[0], noTet[0]);
if (infoRefinement > 1)
{
cout << Form("\n%15s%8d%8d%8d%8d", "removed",
noPt[1]-noPt[0], noEdg[1]-noEdg[0]
, noTr[1]-noTr[0], noTet[1]-noTet[0])
<< Form("\n%15s%8d%8d%8d%8d", "generated",
noPt[2]-noPt[1], noEdg[2]-noEdg[1],
noTr[2]-noTr[1], noTet[2]-noTet[1])
<< Form("\n%15s%8d%8d%8d%8d", "added",
noOfPoints-noPt[2], noOfEdges-noEdg[2],
noOfTriangles-noTr[2], noOfTetrahedra-noTet[2])
<< " for green closure";
}
cout << Form("\n%15s%8d%8d%8d%8d Space : %1.6f MB","total",
noOfPoints, noOfEdges, noOfTriangles,noOfTetrahedra,
float(space)/1e6);
cout << Form("\n%15s%8.2f%8.2f%8.2f%8.2f Blocks: %1d\n","ratio",
float(noOfPoints)/noPt[0], float(noOfEdges)/noEdg[0],
float(noOfTriangles)/noTr[0], float(noOfTetrahedra)/noTet[0],
blocks);
if (infoRefinement > 2)
{
int npDepth, neDepth, ntDepth, step, np=0, ne=0, nt=0, ntet=0;
cout << "\n Grid Statistics:"
<< Form("\n%15s%8s%8s%8s%8s%10s", "depth", "points", "edges",
"tri's","tetra's","(total)");
FORALL(trList, step)
{
npDepth = (step<= ptList.h)? ptList[step]->noOfElements : 0;
neDepth = (step<=edgList.h)? edgList[step]->noOfElements : 0;
ntDepth = (step<= trList.h)? trList[step]->noOfElements : 0;
cout << Form("\n%15d%8d%8d%8d%8d", step,
npDepth, neDepth, ntDepth,
tetList[step]->noOfElements);
np += npDepth;
ne += neDepth;
nt += ntDepth;
ntet +=tetList[step]->noOfElements;
}
cout << Form("\n%15s%8d%8d%8d%8d", "total sum", np, ne, nt, ntet);
cout << "\n";
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void PT3:: print(ostream& os) const
{
os << "Point " << node << " BC " << (int) boundP << " class "
<< (int) classA << " mark " << (int) mark
<< " x: " << x << " y: " << y << " z: " << z << "\n";
}
//-------------------------------------------------------------------------
void PT3:: reset()
{
boundP = INTERIOR;
mark = classA = depth = NUL;
node = 0;
}
//-------------------------------------------------------------------------
void PT3:: getNodeCoord(Vector<Real>& v) const { v[1]=x; v[2]=y; v[3]=z; }
//-------------------------------------------------------------------------
void PT3:: getCoordinates(Vector<Real>& v) const { v[1]=x; v[2]=y; v[3]=z; }
//-------------------------------------------------------------------------
void PT3:: readBC(BufferedParser& parser)
{
readBCValues(&boundP, &classA, parser);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void EDG3:: print(ostream& os) const
{
os << "Edge with points " << p1->node << " " << p2->node
<< " type " << (int) type << " BC " << (int) boundP
<< " class " << (int) classA << " isoP " << (int) isoP
<< " depth " << (int) depth << "\n";
}
void EDG3:: reset()
{
boundP = INTERIOR;
isoP = PLANE;
classA = type = depth = mark = NUL;
p1 = p2 = pm = 0;
node = 0;
father = firstSon = 0;
}
//-------------------------------------------------------------------------
void EDG3:: compNormalVector(Vector<Real>& /*n*/, Real* /*length*/) const
{
notImplemented("EDG3::compNormalVector");
}
//-------------------------------------------------------------------------
// -- point at unit coordinate -> point at x
void EDG3:: 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;
x[2] = L1*p1->y + L2*p2->y;
x[3] = L1*p1->z + L2*p2->z;
}
//-------------------------------------------------------------------------
void EDG3:: getNodeCoord(Vector<Real>& v) const
{
v[1] = 0.5*(p1->x + p2->x);
v[2] = 0.5*(p1->y + p2->y);
v[3] = 0.5*(p1->z + p2->z);
}
//-------------------------------------------------------------------------
void EDG3:: getCoordinates(Vector<Real>& cP1, Vector<Real>& cP2) const
{
p1->getCoordinates(cP1);
p2->getCoordinates(cP2);
}
//-------------------------------------------------------------------------
Real EDG3:: hMax() const
{
return sqrt(sqr(p2->x - p1->x) + sqr(p2->y - p1->y) + sqr(p2->z - p1->z));
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void TR3:: print(ostream& os) const
{
TR::print(os);
os << " type " << (int) type
<< " BC " << (int) boundP
<< " class " << (int) classA << " isoP " << (int) isoP
<< " mark "<< (int)mark << " depth " << (int) depth << "\n";
}
void TR3:: reset()
{
p1 = p2 = p3 = 0;
e1 = e2 = e3 = 0;
mark = classA = depth = NUL;
isoP = PLANE;
boundP = INTERIOR;
type = MESH::T_WHITE;
next = prev = father = firstSon = 0;
node = 0;
t31 = t32 = 0;
}
//-------------------------------------------------------------------------
void TR3:: getEdges(EDG3* eds[]) const
{
eds[1] = e1; eds[2] = e2; eds[3] = e3;
}
//-------------------------------------------------------------------------
// -- point at unit coordinate -> point at x
void TR3:: realCoordinates(const Vector<Real>& u, Vector<Real>& x) const
{
Real L1 = 1.0 - u[1] - u[2];
Real L2 = u[1];
Real L3 = u[2];
x[1] = L1*p1->x + L2*p2->x + L3*p3->x;
x[2] = L1*p1->y + L2*p2->y + L3*p3->y;
x[3] = L1*p1->z + L2*p2->z + L3*p3->z;
}
//-------------------------------------------------------------------------
void TR3:: getCoordinates(Vector<Real>& cP1, Vector<Real>& cP2,
Vector<Real>& cP3) const
{
p1->getCoordinates(cP1);
p2->getCoordinates(cP2);
p3->getCoordinates(cP3);
}
//-------------------------------------------------------------------------
Real TR3:: volume() const
{
Vector<Real> u(3), v(3), w(3);
u[1] = p2->x - p1->x;
u[2] = p2->y - p1->y;
u[3] = p2->z - p1->z;
v[1] = p3->x - p1->x;
v[2] = p3->y - p1->y;
v[3] = p3->z - p1->z;
CrossProduct(w,u,v);
return 0.5 * Norm(w);
}
//-------------------------------------------------------------------------
void TR3:: readBC(BufferedParser& parser)
{
readBCValues(&boundP, &classA, &isoP, parser);
if (isoP > e1->isoP) e1->isoP = isoP;
if (isoP > e2->isoP) e2->isoP = isoP;
if (isoP > e3->isoP) e3->isoP = isoP;
if (boundP == DIRICHLET)
{
e1->boundP = e2->boundP = e3->boundP = boundP;
e1->classA = e2->classA = e3->classA = classA;
p1->boundP = p2->boundP = p3->boundP = boundP;
p1->classA = p2->classA = p3->classA = classA;
}
else if (boundP == CAUCHY || boundP == NEUMANN)
{
if (e1->boundP != DIRICHLET)
{
e1->boundP = boundP;
e1->classA = classA;
}
if (e2->boundP != DIRICHLET)
{
e2->boundP = boundP;
e2->classA = classA;
}
if (e3->boundP != DIRICHLET)
{
e3->boundP = boundP;
e3->classA = classA;
}
if (p1->boundP != DIRICHLET)
{
p1->boundP = boundP;
p1->classA = classA;
}
if (p2->boundP != DIRICHLET)
{
p2->boundP = boundP;
p2->classA = classA;
}
if (p3->boundP != DIRICHLET)
{
p3->boundP = boundP;
p3->classA = classA;
}
}
else
{
cout << "\n*** readBC: Invalid Boundary Type for Face " << boundP << "\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
void TR3:: readHyperBC(BufferedParser& /*parser*/)
{
if (boundP == DIRICHLET)
{
e1->boundP = e2->boundP = e3->boundP = boundP;
e1->classA = e2->classA = e3->classA = classA;
p1->boundP = p2->boundP = p3->boundP = boundP;
p1->classA = p2->classA = p3->classA = classA;
}
else if (boundP == CAUCHY || boundP == NEUMANN)
{
if (e1->boundP != DIRICHLET)
{
e1->boundP = boundP;
e1->classA = classA;
}
if (e2->boundP != DIRICHLET)
{
e2->boundP = boundP;
e2->classA = classA;
}
if (e3->boundP != DIRICHLET)
{
e3->boundP = boundP;
e3->classA = classA;
}
if (p1->boundP != DIRICHLET)
{
p1->boundP = boundP;
p1->classA = classA;
}
if (p2->boundP != DIRICHLET)
{
p2->boundP = boundP;
p2->classA = classA;
}
if (p3->boundP != DIRICHLET)
{
p3->boundP = boundP;
p3->classA = classA;
}
}
else
{
cout << "\n*** readBC: Invalid Boundary Type for Face " << boundP << "\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
void TET3:: compJ(Jacobian& j, Bool /*checkDetJ*/) const
{
j.J(1,1) = p2->x - p1->x; // Jacobian J(i,k) = dx(k)/dr(i)
j.J(2,1) = p3->x - p1->x;
j.J(3,1) = p4->x - p1->x;
j.J(1,2) = p2->y - p1->y;
j.J(2,2) = p3->y - p1->y;
j.J(3,2) = p4->y - p1->y;
j.J(1,3) = p2->z - p1->z;
j.J(2,3) = p3->z - p1->z;
j.J(3,3) = p4->z - p1->z;
j.detJ = detJ;
// = j.J(1,1) * (j.J(2,2)*j.J(3,3) - j.J(3,2)*j.J(2,3))
// + j.J(3,1) * (j.J(1,2)*j.J(2,3) - j.J(2,2)*j.J(1,3))
// + j.J(2,1) * (j.J(1,3)*j.J(3,2) - j.J(3,3)*j.J(1,2));
// if (j.detJ <= 0.0) cout << "\n*** Jacobi-Det. <= 0 " << j.detJ
// << " encountered\n";
}
void TET3:: compJinv(Jacobian& j) const
{
int i, k;
compJ(j);
j.Jinv(1,1) = j.J(2,2)*j.J(3,3) - j.J(3,2)*j.J(2,3);
j.Jinv(1,2) = j.J(1,3)*j.J(3,2) - j.J(1,2)*j.J(3,3);
j.Jinv(1,3) = j.J(1,2)*j.J(2,3) - j.J(1,3)*j.J(2,2);
j.Jinv(2,1) = j.J(2,3)*j.J(3,1) - j.J(3,3)*j.J(2,1);
j.Jinv(2,2) = j.J(1,1)*j.J(3,3) - j.J(3,1)*j.J(1,3);
j.Jinv(2,3) = j.J(1,3)*j.J(2,1) - j.J(2,3)*j.J(1,1);
j.Jinv(3,1) = j.J(2,1)*j.J(3,2) - j.J(3,1)*j.J(2,2);
j.Jinv(3,2) = j.J(1,2)*j.J(3,1) - j.J(3,2)*j.J(1,1);
j.Jinv(3,3) = j.J(1,1)*j.J(2,2) - j.J(2,1)*j.J(1,2);
for (i=1; i<=3; ++i)
for (k=1; k<=3; ++k) j.Jinv(i,k) /= j.detJ;
}
//-------------------------------------------------------------------------
void TET3:: print(ostream& os) const
{
TET::print(os);
os << " type " << (int) type << " class " << (int) classA
<< " mark " << (int)mark << " depth " << (int) depth << "\n";
}
//-------------------------------------------------------------------------
void TET3:: reset()
{
mark = classA = depth = NUL;
boundP = INTERIOR;
detJ = -1.0;
father = firstSon = 0;
node = 0;
type = MESH::T_WHITE;
error = 0.0;
p1 = p2 = p3 = p4 = 0;
n1 = n2 = n3 = n4 = 0;
e1 = e2 = e3 = e4 = e5 = e6 = 0;
}
void TET3:: setBoundary()
{
if ( !(n1->castToTET3() && n2->castToTET3() &&
n3->castToTET3() && n4->castToTET3()) ) boundP = BOUNDARY;
}
//-------------------------------------------------------------------------
void TET3:: getPoints(PT3* pts[]) const
{
pts[1] = p1; pts[2] = p2; pts[3] = p3; pts[4] = p4;
}
//-------------------------------------------------------------------------
void TET3:: getEdges(EDG3* eds[]) const
{
eds[1] = e1; eds[2] = e2; eds[3] = e3;
eds[4] = e4; eds[5] = e5; eds[6] = e6;
}
//-------------------------------------------------------------------------
void TET3:: getEdgeOrientation(Vector<Bool>& orient) const
{
int i;
FORALL(orient,i) orient[i] = False;
if (p1->node == 0 || p2->node == 0 || p3->node == 0 ||p4->node == 0) // !!!
{
cout <<"\n*** TET3:: getEdgeOrientation: node == 0 !\n";
cout.flush(); abort();
}
if (p1->node < p2->node) orient[1] = True;
if (p2->node < p3->node) orient[2] = True;
if (p3->node < p1->node) orient[3] = True;
if (p3->node < p4->node) orient[4] = True;
if (p2->node < p4->node) orient[5] = True;
if (p1->node < p4->node) orient[6] = True;
}
//-------------------------------------------------------------------------
void TET3:: markTriangles(TET3* tet)
{
TET* nbTet;
nbTet = n1->castToTET3();
if ( nbTet == tet) { mark++; return; }
nbTet = n2->castToTET3();
if ( nbTet == tet) { mark += 2; return; }
nbTet = n3->castToTET3();
if ( nbTet == tet) { mark += 4; return; }
nbTet = n4->castToTET3();
if ( nbTet == tet) { mark += 8; return; }
cout << "\n Impossible state in markTriangles\n";
}
//-------------------------------------------------------------------------
void TET3:: getMarkedTriangles(int trsMark[])
{
for (int i=1; i<=4; i++) trsMark[i] = 0;
switch (mark) {
case 0: break;
case 1: trsMark[1] = 1; break;
case 2: trsMark[2] = 1; break;
case 3: trsMark[2] = 1; trsMark[1] = 1; break;
case 4: trsMark[3] = 1; break;
case 5: trsMark[3] = 1; trsMark[1] = 1; break;
case 6: trsMark[3] = 1; trsMark[2] = 1; break;
case 7: trsMark[3] = 1; trsMark[2] = 1; trsMark[1] = 1; break;
case 8: trsMark[4] = 1; break;
case 9: trsMark[4] = 1; trsMark[1] = 1; break;
case 10: trsMark[4] = 1; trsMark[2] = 1; break;
case 11: trsMark[4] = 1; trsMark[2] = 1; trsMark[1] = 1; break;
case 12: trsMark[4] = 1; trsMark[3] = 1; break;
case 13: trsMark[4] = 1; trsMark[3] = 1; trsMark[1] = 1; break;
case 14: trsMark[4] = 1; trsMark[3] = 1; trsMark[2] = 1; break;
case 15: trsMark[4] = 1; trsMark[3] = 1; trsMark[2] = 1;
trsMark[1] = 1; break;
default: cout << "\n dangerous state in getMarkedTriangles\n";
}
return;
}
//-------------------------------------------------------------------------
void TET3:: getAllFaces(Vector<PATCH*>& trs) const
{
TR3 *t1, *t2, *t3, *t4;
TET3 *nt;
void *t;
// t1
if ( (nt = n1->castToTET3()) )
{
t1 = new TR3;
t1->boundP = INTERIOR;
t1->depth = depth;
t1->p1 = p2;
t1->p2 = p3;
t1->p3 = p4;
t1->e1 = e4;
t1->e2 = e5;
t1->e3 = e2;
completeT(t1,0);
// t1->t31 = this; causing warnings, subst. by the next two lines
t = (void*) this;
t1->t31 = (TET3*) t;
t1->t32 = nt;
}
else
{
t1 = n1->castToTR3();
// t1->t31 = this; causing warnings, subst. by the next two lines
t = (void*) this;
t1->t31 = (TET3*) t;
t1->t32 = 0;
}
// t2
if ( (nt = n2->castToTET3()) )
{
t2 = new TR3;
t2->boundP = INTERIOR;
t2->depth = depth;
t2->p1 = p1;
t2->p2 = p4;
t2->p3 = p3;
t2->e1 = e4;
t2->e2 = e1;
t2->e3 = e6;
completeT(t2,0);
// t2->t31 = this; causing warnings, subst. by the next two lines
t = (void*) this;
t2->t31 = (TET3*) t;
t2->t32 = nt;
}
else
{
t2 = n2->castToTR3();
// t2->t31 = this; causing warnings, subst. by the next two lines
t = (void*) this;
t2->t31 = (TET3*) t;
t2->t32 = 0;
}
// t3
if ( (nt = n3->castToTET3()) )
{
t3 = new TR3;
t3->boundP = INTERIOR;
t3->depth = depth;
t3->p1 = p1;
t3->p2 = p2;
t3->p3 = p4;
t3->e1 = e5;
t3->e2 = e6;
t3->e3 = e3;
completeT(t3,0);
// t3->t31 = this; causing warnings, subst. by the next two lines
t = (void*) this;
t3->t31 = (TET3*) t;
t3->t32 = nt;
}
else
{
t3 = n3->castToTR3();
// t3->t31 = this; causing warnings, subst. by the next two lines
t = (void*) this;
t3->t31 = (TET3*) t;
t3->t32 = 0;
}
// t4
if ( (nt = n4->castToTET3()) )
{
t4 = new TR3;
t4->boundP = INTERIOR;
t4->depth = depth;
t4->p1 = p1;
t4->p2 = p3;
t4->p3 = p2;
t4->e1 = e2;
t4->e2 = e3;
t4->e3 = e1;
completeT(t4,0);
// t4->t31 = this; causing warnings, subst. by the next two lines
t = (void*) this;
t4->t31 = (TET3*) t;
t4->t32 = nt;
}
else
{
t4 = n4->castToTR3();
// t4->t31 = this; causing warnings, subst. by the next two lines
t = (void*) this;
t4->t31 = (TET3*) t;
t4->t32 = 0;
}
trs[1] = t1; trs[2] = t2; trs[3] = t3; trs[4] = t4;
}
//-------------------------------------------------------------------------
void TET3:: getTriangles(TR3* trs[])
{
TR3 *t1, *t2, *t3, *t4;
TET3 *nt;
// t1
if ( (nt = n1->castToTET3()) )
{
t1 = new TR3;
t1->boundP = INTERIOR;
t1->depth = depth;
t1->p1 = p2;
t1->p2 = p3;
t1->p3 = p4;
t1->e1 = e4;
t1->e2 = e5;
t1->e3 = e2;
completeT(t1,0);
t1->t31 = this;
t1->t32 = nt;
}
else
{
t1 = n1->castToTR3();
t1->t31 = this;
t1->t32 = 0;
}
// t2
if ( (nt = n2->castToTET3()) )
{
t2 = new TR3;
t2->boundP = INTERIOR;
t2->depth = depth;
t2->p1 = p1;
t2->p2 = p4;
t2->p3 = p3;
t2->e1 = e4;
t2->e2 = e1;
t2->e3 = e6;
completeT(t2,0);
t2->t31 = this;
t2->t32 = nt;
}
else
{
t2 = n2->castToTR3();
t2->t31 = this;
t2->t32 = 0;
}
// t3
if ( (nt = n3->castToTET3()) )
{
t3 = new TR3;
t3->boundP = INTERIOR;
t3->depth = depth;
t3->p1 = p1;
t3->p2 = p2;
t3->p3 = p4;
t3->e1 = e5;
t3->e2 = e6;
t3->e3 = e3;
completeT(t3,0);
t3->t31 = this;
t3->t32 = nt;
}
else
{
t3 = n3->castToTR3();
t3->t31 = this;
t3->t32 = 0;
}
// t4
if ( (nt = n4->castToTET3()) )
{
t4 = new TR3;
t4->boundP = INTERIOR;
t4->depth = depth;
t4->p1 = p1;
t4->p2 = p3;
t4->p3 = p2;
t4->e1 = e2;
t4->e2 = e3;
t4->e3 = e1;
completeT(t4,0);
t4->t31 = this;
t4->t32 = nt;
}
else
{
t4 = n4->castToTR3();
t4->t31 = this;
t4->t32 = 0;
}
trs[1] = t1; trs[2] = t2; trs[3] = t3; trs[4] = t4;
}
//-------------------------------------------------------------------------
void TET3:: returnFaces(Vector<PATCH*>& trs) const
{
for (int i=1; i<=4; i++) if (!(trs[i]->onBoundary())) delete trs[i];
}
//-------------------------------------------------------------------------
void TET3:: returnTriangles(TR3* trs[]) const
{
for (int i=1; i<=4; i++) if (!(trs[i]->onBoundary())) delete trs[i];
}
//-------------------------------------------------------------------------
void TET3:: getNeighbours(Vector<PATCH*>& neighbours) const
{
neighbours[1] = n1->castToTET3();
neighbours[2] = n2->castToTET3();
neighbours[3] = n3->castToTET3();
neighbours[4] = n4->castToTET3();
}
//-------------------------------------------------------------------------
// -- point at unit coordinate -> point at x
void TET3:: realCoordinates(const Vector<Real>& u, Vector<Real>& x) const
{
Real L1 = 1.0 - u[1] - u[2] - u[3];
Real L2 = u[1];
Real L3 = u[2];
Real L4 = u[3];
x[1] = L1*p1->x + L2*p2->x + L3*p3->x + L4*p4->x;
x[2] = L1*p1->y + L2*p2->y + L3*p3->y + L4*p4->y;
x[3] = L1*p1->z + L2*p2->z + L3*p3->z + L4*p4->z;
}
//-------------------------------------------------------------------------
// -- point at x -> unit coordinates
void TET3:: unitCoordinates(const Vector<Real>& x, Vector<Real>& u) const
{
int i, k;
Jacobian J;
Vector<Real> xn(3);
compJinv(J);
xn[1] = x[1] - p1->x;
xn[2] = x[2] - p1->y;
xn[3] = x[3] - p1->z;
for (i=1; i<=3; ++i)
{
u[i] = 0.0;
for (k=1; k<=3; ++k) u[i] += J.Jinv(k,i) * xn[k];
} // Jacobian is transposed!
}
//-------------------------------------------------------------------------
void TET3:: getCoordinates(Vector<Real>& cP1, Vector<Real>& cP2,
Vector<Real>& cP3, Vector<Real>& cP4) const
{
p1->getCoordinates(cP1);
p2->getCoordinates(cP2);
p3->getCoordinates(cP3);
p4->getCoordinates(cP4);
}
//-------------------------------------------------------------------------
Real TET3:: volume() const { return oneSixth*Abs(detJ); }
//-------------------------------------------------------------------------
int TET3:: NoOfSons() const
{
if (type == MESH::T_RED) return 8;
else if (type == MESH::T_GREEN_1) return 2;
else if (type == MESH::T_GREEN_2A) return 4;
else if (type == MESH::T_GREEN_2B) return 3;
else if (type == MESH::T_GREEN_3) return 4;
else
{
cout << "\n*** TET3::NoOfSons(): unknown type " << type << "\n";
cout.flush(); abort();
return 0;
}
}
//-------------------------------------------------------------------------
Bool TET3:: isGreen() const
{
switch(type)
{
case MESH::T_GREEN:
case MESH::T_GREEN_1:
case MESH::T_GREEN_2A:
case MESH::T_GREEN_2B:
case MESH::T_GREEN_3: return True;
default: return False;
}
}
Bool TR3:: isGreen() const
{
switch(type)
{
case MESH::T_GREEN:
case MESH::T_GREEN_1:
case MESH::T_GREEN_2A:
case MESH::T_GREEN_2B:
case MESH::T_GREEN_3: return True;
default: return False;
}
}
Bool EDG3:: isGreen() const
{
switch(type)
{
case MESH::T_GREEN:
case MESH::T_GREEN_1:
case MESH::T_GREEN_2A:
case MESH::T_GREEN_2B:
case MESH::T_GREEN_3: return True;
default: return False;
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void EFTET3:: reset()
{
for (int i=0; i<4; ++i) faceNode[i] = 0;
TET3::reset();
}
void EFTET3:: getFaceNodes(Vector<int>& nodes) const
{
int i;
for (i=1; i<=4; ++i) nodes[i] = faceNode[i-1];
}
int& EFTET3:: getFaceNode(int face) { return faceNode[face-1]; }
int& EFTET3:: getMyFaceNode(PATCH* patch)
{
if (n1 == patch) return faceNode[0];
if (n2 == patch) return faceNode[1];
if (n3 == patch) return faceNode[2];
if (n4 == patch) return faceNode[3];
cout << "\n*** EFTET3:: getMyFaceNode: patch not found\n";
cout.flush(); abort();
return faceNode[0]; // dummy return
}
//-------------------------------------------------------------------------
/*-------------------------------------------------------------------------*/
void MESH3:: setXYZpmCircle(PT3* pm, const EDG3* ed)
{
Real midPointX, midPointY, midPointZ, radius_new;
Real hx1, hy1, hz1, hx2, hy2, hz2, radius;
PT3 *P1 = ed->p1, *P2 = ed->p2, *p = pm;
// P1,P2 on surface of the ball
midPointX = 0.0;
midPointY = 0.0;
midPointZ = 0.0;
hx1 = P1->x-midPointX; hy1 = P1->y-midPointY; hz1 = P1->z-midPointZ;
radius = sqrt(hx1*hx1+hy1*hy1+hz1*hz1);
p->x = (P1->x)+((P2->x)-(P1->x))*HALF;
p->y = (P1->y)+((P2->y)-(P1->y))*HALF;
p->z = (P1->z)+((P2->z)-(P1->z))*HALF;
hx2 = p->x-midPointX; hy2 = p->y-midPointY; hz2 = p->z-midPointZ;
radius_new = sqrt(hx2*hx2+hy2*hy2+hz2*hz2);
p->x = midPointX + radius*hx2/radius_new;
p->y = midPointY + radius*hy2/radius_new;
p->z = midPointZ + radius*hz2/radius_new;
}
/*-------------------------------------------------------------------------*/
void MESH3:: setXYZpmCyl(PT3* pm, const EDG3* ed)
{
Real radius,hx1,hy1,hx2,hy2;
Real radius_new;
PT3 *P1 = ed->p1, *P2 = ed->p2;
hx1 = P1->x; hy1 = P1->y;
radius = sqrt(hx1*hx1+hy1*hy1);
pm->x = (P1->x)+((P2->x)-(P1->x))*HALF;
pm->y = (P1->y)+((P2->y)-(P1->y))*HALF;
pm->z = (P1->z)+((P2->z)-(P1->z))*HALF;
hx2 = pm->x; hy2 = pm->y;
radius_new = sqrt(hx2*hx2+hy2*hy2);
pm->x = radius*hx2/radius_new;
pm->y = radius*hy2/radius_new;
}
//-------------------------------------------------------------------------
void MESH3:: resetElemIter(Bool lastStep)
{
if (elemIter != 0)
{
cout << "\n*** resetElemIter: elemIter != 0"
<< " \t (previous operation not completed?)\n";
cout.flush(); abort();
}
if (lastStep) iterLevel = tetList.high();
else iterLevel = 0;
}
//-------------------------------------------------------------------------
PATCH* MESH3:: elemIterAll()
{
if (elemIter != 0) // iterate rest of list
{
for (elemIter=elemIter->Next(); elemIter; elemIter=elemIter->Next())
if (!elemIter->refined()) return elemIter;
if (++iterLevel > tetList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (elemIter=tetList[iterLevel]->first; elemIter;
elemIter=elemIter->Next())
if (!elemIter->refined()) return elemIter;
if (++iterLevel > tetList.high()) return 0;
}
}
//-------------------------------------------------------------------------
PATCH* MESH3:: elemIterAllOfHistory()
{
if (elemIter != 0) // iterate rest of list
{
for (elemIter=elemIter->Next(); elemIter; elemIter=elemIter->Next())
return elemIter;
if (++iterLevel > tetList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (elemIter=tetList[iterLevel]->first; elemIter;
elemIter=elemIter->Next()) return elemIter;
if (++iterLevel > tetList.high()) return 0;
}
}
//-------------------------------------------------------------------------
PATCH* MESH3:: 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* MESH3:: 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* MESH3:: 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