HPREF_ELEMENT_TYPE ClassifyTet(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
{
int ep1(0), ep2(0), ep3(0), ep4(0), cp1(0), cp2(0), cp3(0), cp4(0), fp1, fp2, fp3, fp4;
int isedge1(0), isedge2(0), isedge3(0), isedge4(0), isedge5(0), isedge6(0);
int isfedge1, isfedge2, isfedge3, isfedge4, isfedge5, isfedge6;
bool isface[4];
HPREF_ELEMENT_TYPE type = HP_NONE;
int debug = 0;
for (int j = 0;j < 4; j++)
{
if (el.pnums[j] == 444) debug++;
if (el.pnums[j] == 115) debug++;
if (el.pnums[j] == 382) debug++;
if (el.pnums[j] == 281) debug++;
}
if (debug < 4) debug = 0;
// *testout << "new el" << endl;
for (int j = 0; j < 4; j++)
for (int k = 0; k < 4; k++)
{
if (j == k) continue;
if (type) break;
int pi3 = 0;
while (pi3 == j || pi3 == k) pi3++;
int pi4 = 6 - j - k - pi3;
// preserve orientation
int sort[4];
sort[0] = j; sort[1] = k; sort[2] = pi3; sort[3] = pi4;
int cnt = 0;
for (int jj = 0; jj < 4; jj++)
for (int kk = 0; kk < 3; kk++)
if (sort[kk] > sort[kk+1])
{
cnt++;
Swap (sort[kk], sort[kk+1]);
}
if (cnt % 2 == 1) Swap (pi3, pi4);
ep1 = edgepoint.Test (el.pnums[j]);
ep2 = edgepoint.Test (el.pnums[k]);
ep3 = edgepoint.Test (el.pnums[pi3]);
ep4 = edgepoint.Test (el.pnums[pi4]);
cp1 = cornerpoint.Test (el.pnums[j]);
cp2 = cornerpoint.Test (el.pnums[k]);
cp3 = cornerpoint.Test (el.pnums[pi3]);
cp4 = cornerpoint.Test (el.pnums[pi4]);
isedge1 = edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[k]));
isedge2 = edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[pi3]));
isedge3 = edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[pi4]));
isedge4 = edges.Used (INDEX_2::Sort (el.pnums[k], el.pnums[pi3]));
isedge5 = edges.Used (INDEX_2::Sort (el.pnums[k], el.pnums[pi4]));
isedge6 = edges.Used (INDEX_2::Sort (el.pnums[pi3], el.pnums[pi4]));
if (debug)
{
cout << "debug" << endl;
*testout << "debug" << endl;
*testout << "ep = " << ep1 << ep2 << ep3 << ep4 << endl;
*testout << "cp = " << cp1 << cp2 << cp3 << cp4 << endl;
*testout << "edge = " << isedge1 << isedge2 << isedge3 << isedge4 << isedge5 << isedge6 << endl;
}
for (int j = 0; j < 4; j++) isface[j] = false;
for (int l = 0; l < 4; l++)
{
INDEX_3 i3(0,0,0);
switch (l)
{
case 0: i3.I1() = el.pnums[k]; i3.I2() = el.pnums[pi3]; i3.I3() = el.pnums[pi4]; break;
case 1: i3.I1() = el.pnums[j]; i3.I2() = el.pnums[pi3]; i3.I3() = el.pnums[pi4]; break;
case 2: i3.I1() = el.pnums[j]; i3.I2() = el.pnums[k]; i3.I3() = el.pnums[pi4]; break;
case 3: i3.I1() = el.pnums[j]; i3.I2() = el.pnums[k]; i3.I3() = el.pnums[pi3]; break;
}
i3.Sort();
if (faces.Used (i3))
{
int domnr = faces.Get(i3);
if (domnr == -1 || domnr == el.GetIndex())
isface[l] = true;
}
}
/*
isface1 = faces.Used (INDEX_3::Sort (el.pnums[k], el.pnums[pi3], el.pnums[pi4]));
isface2 = faces.Used (INDEX_3::Sort (el.pnums[j], el.pnums[pi3], el.pnums[pi4]));
isface3 = faces.Used (INDEX_3::Sort (el.pnums[j], el.pnums[k], el.pnums[pi4]));
isface4 = faces.Used (INDEX_3::Sort (el.pnums[j], el.pnums[k], el.pnums[pi3]));
*/
isfedge1 = isfedge2 = isfedge3 = isfedge4 = isfedge5 = isfedge6 = 0;
for (int l = 0; l < 6; l++)
{
INDEX_2 i2(0,0);
switch (l)
{
case 0: i2.I1() = el.pnums[j]; i2.I2() = el[k]; break;
case 1: i2.I1() = el.pnums[j]; i2.I2() = el.pnums[pi3]; break;
case 2: i2.I1() = el.pnums[j]; i2.I2() = el.pnums[pi4]; break;
case 3: i2.I1() = el.pnums[k]; i2.I2() = el.pnums[pi3]; break;
case 4: i2.I1() = el.pnums[k]; i2.I2() = el.pnums[pi4]; break;
case 5: i2.I1() = el.pnums[pi3]; i2.I2() = el.pnums[pi4]; break;
}
i2.Sort();
if (face_edges.Used (i2))
{
int domnr = face_edges.Get(i2);
if (domnr == -1 || domnr == el.GetIndex())
{
switch (l)
{
case 0: isfedge1 = 1; break;
case 1: isfedge2 = 1; break;
case 2: isfedge3 = 1; break;
case 3: isfedge4 = 1; break;
case 4: isfedge5 = 1; break;
case 5: isfedge6 = 1; break;
}
}
}
}
/*
isfedge1 = face_edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[k]));
isfedge2 = face_edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[pi3]));
isfedge3 = face_edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[pi4]));
isfedge4 = face_edges.Used (INDEX_2::Sort (el.pnums[k], el.pnums[pi3]));
isfedge5 = face_edges.Used (INDEX_2::Sort (el.pnums[k], el.pnums[pi4]));
isfedge6 = face_edges.Used (INDEX_2::Sort (el.pnums[pi3], el.pnums[pi4]));
*/
fp1 = fp2 = fp3 = fp4 = 0;
for (int l = 0; l < 4; l++)
{
int pti(0);
switch (l)
{
case 0: pti = el.pnums[j]; break;
case 1: pti = el.pnums[k]; break;
case 2: pti = el.pnums[pi3]; break;
case 3: pti = el.pnums[pi4]; break;
}
int domnr = facepoint[pti];
if (domnr == -1 || domnr == el.GetIndex())
{
switch (l)
{
case 0: fp1 = 1; break;
case 1: fp2 = 1; break;
case 2: fp3 = 1; break;
case 3: fp4 = 1; break;
}
}
}
/*
ep1 |= cp1;
ep2 |= cp2;
ep3 |= cp3;
ep4 |= cp4;
fp1 |= ep1;
fp2 |= ep2;
fp3 |= ep3;
fp4 |= ep4;
*/
/*
fp1 = facepoint[el.pnums[j]] != 0;
fp2 = facepoint[el.pnums[k]] != 0;
fp3 = facepoint[el.pnums[pi3]] != 0;
fp4 = facepoint[el.pnums[pi4]] != 0;
*/
// cout << "marked faces: "
// << isface[0] << isface[1] << isface[2] << isface[3]
// << ", num = " << isface[0]+isface[1]+isface[2]+isface[3] << endl;
bool sp1 = cp1
|| (ep1 && !isedge1 && !isedge2 && !isedge3)
|| (fp1 && !isfedge1 && !isfedge2 && !isfedge3);
bool sp2 = cp2
|| (ep2 && !isedge1 && !isedge4 && !isedge5)
|| (fp2 && !isfedge1 && !isfedge4 && !isfedge5);
bool sp3 = cp3
|| (ep3 && !isedge2 && !isedge4 && !isedge6)
|| (fp3 && !isfedge2 && !isfedge4 && !isfedge6);
bool sp4 = cp4
|| (ep4 && !isedge3 && !isedge5 && !isedge6)
|| (fp4 && !isfedge3 && !isfedge5 && !isfedge6);
bool se1 = isedge1 || (isfedge1 && !isface[2] && !isface[3]);
bool se2 = isedge2 || (isfedge2 && !isface[1] && !isface[3]);
bool se3 = isedge3 || (isfedge3 && !isface[1] && !isface[2]);
bool se4 = isedge4 || (isfedge4 && !isface[0] && !isface[3]);
bool se5 = isedge5 || (isfedge5 && !isface[0] && !isface[2]);
bool se6 = isedge6 || (isfedge6 && !isface[0] && !isface[1]);
// *testout << "sp = " << sp1 << sp2 << sp3 << sp4 << endl;
// *testout << "se = " << se1 << se2 << se3 << se4 << se5 << se6 << endl;
// *testout << "sf = " << isface[0] << isface[1] << isface[2] << isface[3] << endl;
switch (isface[0]+isface[1]+isface[2]+isface[3])
{
case 0:
{
isedge1 |= isfedge1;
isedge2 |= isfedge2;
isedge3 |= isfedge3;
isedge4 |= isfedge4;
isedge5 |= isfedge5;
isedge6 |= isfedge6;
ep1 |= fp1;
ep2 |= fp2;
ep3 |= fp3;
ep4 |= fp4;
switch (isedge1+isedge2+isedge3+isedge4+isedge5+isedge6)
{
case 0:
{
if (!sp1 && !sp2 && !sp3 && !sp4)
type = HP_TET;
if (sp1 && !sp2 && !sp3 && !sp4)
type = HP_TET_0E_1V;
if (sp1 && sp2 && !sp3 && !sp4)
type = HP_TET_0E_2V;
if (sp1 && sp2 && sp3 && !sp4)
type = HP_TET_0E_3V;
if (sp1 && sp2 && sp3 && sp4)
type = HP_TET_0E_4V;
break;
}
case 1:
{
if (!isedge1) break;
if (!sp1 && !sp2 && !sp3 && !sp4)
type = HP_TET_1E_0V;
if (sp1 && !sp2 && !sp3 && !sp4)
type = HP_TET_1E_1VA;
if (!sp1 && !sp2 && !sp3 && sp4)
type = HP_TET_1E_1VB;
if (sp1 && sp2 && !sp3 && !sp4)
type = HP_TET_1E_2VA;
if (sp1 && !sp2 && sp3 && !sp4)
type = HP_TET_1E_2VB;
if (sp1 && !sp2 && !sp3 && sp4)
type = HP_TET_1E_2VC;
if (!sp1 && !sp2 && sp3 && sp4)
type = HP_TET_1E_2VD;
if (sp1 && sp2 && sp3 && !sp4)
type = HP_TET_1E_3VA;
if (sp1 && !sp2 && sp3 && sp4)
type = HP_TET_1E_3VB;
if (sp1 && sp2 && sp3 && sp4)
type = HP_TET_1E_4V;
break;
}
case 2:
{
if (isedge1 && isedge2)
{
if (!sp2 && !sp3 && !sp4)
type = HP_TET_2EA_0V;
if (sp2 && !sp3 && !sp4)
type = HP_TET_2EA_1VA;
if (!sp2 && sp3 && !sp4)
type = HP_TET_2EA_1VB;
if (!sp2 && !sp3 && sp4)
type = HP_TET_2EA_1VC;
if (sp2 && sp3 && !sp4)
type = HP_TET_2EA_2VA;
if (sp2 && !sp3 && sp4)
type = HP_TET_2EA_2VB;
if (!sp2 && sp3 && sp4)
type = HP_TET_2EA_2VC;
if (sp2 && sp3 && sp4)
type = HP_TET_2EA_3V;
}
if (isedge1 && isedge6)
{
if (!sp1 && !sp2 && !sp3 && !sp4)
type = HP_TET_2EB_0V;
if (sp1 && !sp2 && !sp3 && !sp4)
type = HP_TET_2EB_1V;
if (sp1 && sp2 && !sp3 && !sp4)
type = HP_TET_2EB_2VA;
if (sp1 && !sp2 && sp3 && !sp4)
type = HP_TET_2EB_2VB;
if (sp1 && !sp2 && !sp3 && sp4)
type = HP_TET_2EB_2VC;
if (sp1 && sp2 && sp3 && !sp4)
type = HP_TET_2EB_3V;
if (sp1 && sp2 && sp3 && sp4)
type = HP_TET_2EB_4V;
}
break;
}
case 3:
{
if (isedge1 && isedge2 && isedge3)
{
if (!sp2 && !sp3 && !sp4)
type = HP_TET_3EA_0V;
if (sp2 && !sp3 && !sp4)
type = HP_TET_3EA_1V;
if (sp2 && sp3 && !sp4)
type = HP_TET_3EA_2V;
if (sp2 && sp3 && sp4)
type = HP_TET_3EA_3V;
}
if (isedge1 && isedge3 && isedge4)
{
if (!sp3 && !sp4)
type = HP_TET_3EB_0V;
if (sp3 && !sp4)
type = HP_TET_3EB_1V;
if (sp3 && sp4)
type = HP_TET_3EB_2V;
}
if (isedge1 && isedge2 && isedge5)
{
if (!sp3 && !sp4)
type = HP_TET_3EC_0V;
if (sp3 && !sp4)
type = HP_TET_3EC_1V;
if (sp3 && sp4)
type = HP_TET_3EC_2V;
}
if (isedge1 && isedge2 && isedge4)
{
if (!sp4)
type = HP_TET_3ED_3V; // a loop
}
break;
}
}
break;
}
case 1: // one singular face
{
if (!isface[0]) break;
/*
cout << "1F and 1E, isedge = " << isedge1 << isedge2 << isedge3 << isedge4 << isedge5 << isedge6 << endl;
cout << "spoints = " << sp1 << sp2 << sp3 << sp4 << endl;
cout << "cpoints = " << cp1 << cp2 << cp3 << cp4 << endl;
cout << "epoints = " << ep1 << ep2 << ep3 << ep4 << endl;
cout << "fpoints = " << fp1 << fp2 << fp3 << fp4 << endl;
*/
isedge1 |= isfedge1;
isedge2 |= isfedge2;
isedge3 |= isfedge3;
// switch (isedge1+isedge2+isedge3+isedge4+isedge5+isedge6)
switch (se1+se2+se3+se4+se5+se6)
{
case 0:
{
if (!fp1 && !ep2 && !ep3 && !ep4)
type = HP_TET_1F_0E_0V;
if (fp1 && !ep2 && !ep3 && !ep4)
type = HP_TET_1F_0E_1VB;
if (!fp1 && ep2 && !ep3 & !ep4)
type = HP_TET_1F_0E_1VA;
if (!fp1 && ep2 && ep3 & !ep4)
type = HP_TET_1F_0E_2V;
if (!sp1 && sp2 && sp3 && sp4)
type = HP_TET_1F_0E_3V;
break;
}
case 1:
{
if (se1)
{
if (!sp1 && !sp3 && !sp4)
type = HP_TET_1F_1EA_0V;
if (!sp1 && sp2 && sp3 && !sp4)
type = HP_TET_1F_1E_2VA;
if (!sp1 && sp2 && !sp3 && sp4)
type = HP_TET_1F_1E_2VB;
if (!sp1 && !sp2 && sp3 && sp4)
type = HP_TET_1F_1E_2VC;
if (!sp1 && sp2 && sp3 && sp4)
type = HP_TET_1F_1EA_3V;
}
if (se4) // V2-V3
{
if (!sp1 && !sp2 && !sp3 && !sp4)
type = HP_TET_1F_1EB_0V;
if (!sp1 && sp2 && !sp3 && !sp4)
type = HP_TET_1F_1E_1VA;
if (!sp1 && sp2 && sp3 && sp4)
type = HP_TET_1F_1E_3V;
}
if (se5) // V2-V4
{
if (!sp1 && sp2 && !sp3 && !sp4)
type = HP_TET_1F_1E_1VB;
}
break;
}
case 2:
{
if (isedge1 && isedge2)
{
if (sp1 && sp2 && sp3 && !sp4)
type = HP_TET_1F_2Eoo_3V;
}
if (isedge6 && isedge3)
if (!cp1 && !cp2 && !cp3)
type = HP_TET_1F_2E_0VA;
if (isedge6 && isedge2)
{
if (!cp1 && !cp2 && !cp4)
type = HP_TET_1F_2E_0VB;
}
if (se4 && se5)
{ // 2 edges in face
if (!sp1 && sp2 && !sp3 && !sp4)
type = HP_TET_1F_2E_1V;
if (!sp1 && sp2 && sp3 && sp4)
type = HP_TET_1F_2E_3V;
}
break;
}
default:
;
}
break;
}
case 2: // two singular faces
{
if (!isface[0] || !isface[1]) break;
switch (isfedge1+isedge2+isedge3+isedge4+isedge5)
{
case 0:
{
if (!ep1 && !ep2 && !cp3 && !cp4)
{
type = HP_TET_2F_0E_0V;
break;
}
if (!ep1 && !ep2 && !cp3 && cp4)
{
type = HP_TET_2F_0E_1V;
break;
}
break;
}
case 1:
{
// *testout << "so far: 2F, 1E, sp = " << sp1 << sp2 << sp3 << sp4 << endl;
if (isedge4)
{
if (!ep1 && !cp2 && !cp4)
{
type = HP_TET_2F_1E_0VA;
break;
}
if (!sp1 && sp2 && sp3 && sp4)
{
type = HP_TET_2F_1E_3VA;
break;
}
if (sp1 && sp2 && sp3 && sp4)
{
type = HP_TET_2F_1E_4VA;
break;
}
}
if (isedge5 && !ep1 && !cp2 && !cp3)
{
type = HP_TET_2F_1E_0VB;
break;
}
break;
}
default:
*testout << "2F, 2E or more not implemented so far" << endl;
}
break;
}
case 3:
{
if (!isface[3])
if (!cp1 && !cp2 && !cp3)
{
type = HP_TET_3F_0E_0V;
break;
}
break;
}
case 4:
{
*testout << "4 singular faces" << endl;
}
}
if (type != HP_NONE)
{
int pnums[4];
pnums[0] = el.pnums[j];
pnums[1] = el.pnums[k];
pnums[2] = el.pnums[pi3];
pnums[3] = el.pnums[pi4];
for(k=0;k<4;k++) el.pnums[k] = pnums[k];
break;
}
}
if (debug) cout << "type = " << type << endl;
if (type == HP_NONE)
{
// cnt_undef++;
(*testout) << "unclassified element "
<< el.pnums[0] << " "
<< el.pnums[1] << " "
<< el.pnums[2] << " "
<< el.pnums[3] << endl
<< "cp = " << cp1 << cp2 << cp3 << cp4 << endl
<< "ep = " << ep1 << ep2 << ep3 << ep4 << endl
<< "fp = " << fp1 << fp2 << fp3 << fp4 << endl
<< "isedge = " << isedge1 << isedge2 << isedge3
<< isedge4 << isedge5 << isedge6 << endl
<< "isfedge = " << isfedge1 << isfedge2 << isfedge3
<< isfedge4 << isfedge5 << isfedge6 << endl
<< "isface = " << isface[0] << isface[1] << isface[2] << isface[3] << endl;
cout << "unclassified element !!! " << endl;
}
return(type);
}
HPREF_ELEMENT_TYPE ClassifyPrism(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
{
HPREF_ELEMENT_TYPE type = HP_NONE;
int p[6];
for(int m=1;m<=6;m++)
{
int point_sing[6]={0,0,0,0,0,0};
int face_sing[5]={0,0,0,0,0};
int edge_sing[9]={0,0,0,0,0,0,0,0,0};
if(m<4)
{
p[0]= m; p[1]=m%3+1; p[2]=(m%3+1)%3+1;
for(int l=3;l<6;l++) p[l]=p[l-3]+3;
}
else
{
p[0] = m; p[1]=(m%3+1)%3+4; p[2]=m%3+4;
for(int l=3;l<6;l++) p[l]=p[l-3]-3;
}
for(int j=0;j<6;j++)
{
if(cornerpoint.Test(el.PNum(p[j]))) { point_sing[p[j]-1]=3;}
else if(edgepoint.Test(el.PNum(p[j]))) point_sing[p[j]-1]=2;
else if (facepoint[el.PNum(p[j])] == -1 || facepoint[el.PNum(p[j])] == el.GetIndex())
point_sing[p[j]-1] = 1;
}
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (PRISM);
for(int k=0;k<9;k++)
{
INDEX_2 i2 = INDEX_2 :: Sort(el.PNum(p[eledges[k][0]-1]),el.PNum(p[eledges[k][1]-1]));
if (edges.Used(i2)) edge_sing[k] = 2;
else edge_sing[k] = face_edges.Used(i2);
}
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (PRISM);
for (int k=0;k<5;k++)
{
INDEX_3 i3;
if(k<2)
i3 = INDEX_3::Sort(el.pnums[p[elfaces[k][0]-1]-1], el.pnums[p[elfaces[k][1]-1]-1],
el.pnums[p[elfaces[k][2]-1]-1]);
else
{
INDEX_4 i4 = INDEX_4(el.pnums[p[elfaces[k][0]-1]-1], el.pnums[p[elfaces[k][1]-1]-1], el.pnums[p[elfaces[k][2]-1]-1],el.pnums[p[elfaces[k][3]-1]-1]);
i4.Sort();
i3 = INDEX_3(i4.I1(), i4.I2(), i4.I3());
}
if (faces.Used (i3))
{
int domnr = faces.Get(i3);
if (domnr == -1 || domnr == el.GetIndex())
face_sing[k] = 1;
}
}
if (face_sing[1] > face_sing[0]) {m=m+2; continue;}
//int cp = 0;
int qfsing = face_sing[2] + face_sing[3] + face_sing[4];
int tfsing = face_sing[0] + face_sing[1];
int evsing = edge_sing[6] + edge_sing[7] + edge_sing[8];
int ehsing = edge_sing[0] + edge_sing[1] + edge_sing[2] + edge_sing[3] + edge_sing[4] + edge_sing[5];
if (qfsing + tfsing + evsing + ehsing == 0)
{ type = HP_PRISM; break;}
HPREF_ELEMENT_TYPE types[] = {HP_NONE,HP_NONE,HP_NONE};
int fb = (1-face_sing[4])* face_sing[3] * (face_sing[2] + face_sing[3]) + 3*face_sing[4]*face_sing[3]*face_sing[2];
int sve[3] = {edge_sing[7] , edge_sing[8], edge_sing[6]};
if(fb!=qfsing) continue;
switch(fb)
{
case 0:
if (evsing == 0 && ehsing==3*tfsing)
{
types[0] = HP_PRISM;
types[1] = HP_PRISM_1FA_0E_0V;
types[2] = HP_PRISM_2FA_0E_0V;
}
if(evsing > 0 && sve[0] == evsing) // 1 vertical edge 1-4
{
types[0] = HP_PRISM_SINGEDGE;
types[1] = HP_PRISM_1FA_1E_0V;
types[2] = HP_PRISM_2FA_1E_0V;
}
if(sve[0] > 0 && sve[1] > 0 && sve[2] == 0)
{
types[0] = HP_PRISM_SINGEDGE_V12;
types[1] = HP_PRISM_1FA_2E_0V;
types[2] = HP_PRISM_2FA_2E_0V;
}
if(sve[0] > 0 && sve[1] > 0 && sve[2] > 0)
{
types[0] = HP_PRISM_3E_0V;
types[1] = HP_PRISM_1FA_3E_0V;
types[2] = HP_PRISM_2FA_3E_0V;
if ( edge_sing[0] > 1 && edge_sing[2] > 1 &&
edge_sing[4] > 1 && edge_sing[5] > 1 && tfsing==0)
types[0] = HP_PRISM_3E_4EH;
}
break;
case 1:
if(sve[0] <= 1 && sve[1] <= 1)
{
if(sve[2]==0)
{
types[0] = HP_PRISM_1FB_0E_0V;
types[1] = HP_PRISM_1FA_1FB_0E_0V;
types[2] = HP_PRISM_2FA_1FB_0E_0V;
}
else
{
types[0] = HP_PRISM_1FB_1EC_0V;
types[1] = HP_PRISM_1FA_1FB_1EC_0V;
types[2] = HP_PRISM_2FA_1FB_1EC_0V;
}
}
if(sve[0] > 1 && sve[2] >= 1 && sve[1] <= 1)
{
types[0] = HP_PRISM_1FB_2EB_0V;
types[1] = HP_PRISM_1FA_1FB_2EB_0V;
types[2] = HP_PRISM_2FA_1FB_2EB_0V;
}
if(sve[0] > 1 && sve[1] <= 1 && sve[2] == 0) // ea && !eb
{
types[0] = HP_PRISM_1FB_1EA_0V;
types[1] = HP_PRISM_1FA_1FB_1EA_0V;
types[2] = HP_PRISM_2FA_1FB_1EA_0V;
}
if(sve[0] <= 1 && sve[1] > 1 && sve[2] == 0)
types[1] = HP_PRISM_1FA_1FB_1EB_0V;
if(sve[0] > 1 && sve[1]>1)
if(sve[2] == 0) // ea && eb
{
types[0] = HP_PRISM_1FB_2EA_0V;
types[1] = HP_PRISM_1FA_1FB_2EA_0V;
types[2] = HP_PRISM_2FA_1FB_2EA_0V;
}
if(sve[0] <= 1 && sve[1] > 1 && sve[2] >0)
types[1] = HP_PRISM_1FA_1FB_2EC_0V;
if(sve[0] > 1 && sve[1] > 1 && sve[2] >= 1) //sve[2] can also be a face-edge
{
types[0] = HP_PRISM_1FB_3E_0V;
types[1] = HP_PRISM_1FA_1FB_3E_0V;
types[2] = HP_PRISM_2FA_1FB_3E_0V;
}
break;
case 2:
if(sve[0] <= 1)
cout << " **** WARNING: Edge between to different singular faces should be marked singular " << endl;
if(sve[1] <= 1)
if(sve[2] <=1)
{
types[0] = HP_PRISM_2FB_0E_0V;
types[1] = HP_PRISM_1FA_2FB_0E_0V;
types[2] = HP_PRISM_2FA_2FB_0E_0V;
}
else
{
types[0] = HP_PRISM_2FB_1EC_0V;
types[1] = HP_PRISM_1FA_2FB_1EC_0V;
types[2] = HP_PRISM_2FA_2FB_1EC_0V;
}
else
if(sve[2] <= 1)
types[1] = HP_PRISM_1FA_2FB_1EB_0V;
else
{
types[0] = HP_PRISM_2FB_3E_0V;
types[1] = HP_PRISM_1FA_2FB_3E_0V;
types[2] = HP_PRISM_2FA_2FB_3E_0V;
}
break;
case 3:
types[0] = HP_PRISM_3FB_0V;
types[1] = HP_PRISM_1FA_3FB_0V;
types[2] = HP_PRISM_2FA_3FB_0V;
break;
}
type = types[tfsing];
if(type != HP_NONE)
break;
}
/*
*testout << " Prism with pnums " << endl;
for(int j=0;j<6;j++) *testout << el.pnums[j] << "\t";
*testout << endl;
*/
if(type != HP_NONE)
{
int pnums[6];
for(int j=0;j<6;j++) pnums[j] = el.PNum (p[j]);
for(int k=0;k<6;k++) el.pnums[k] = pnums[k];
}
/* *testout << " Classified Prism with pnums " << endl;
for(int j=0;j<6;j++) *testout << el.pnums[j] << "\t";
*testout << endl;
*/
return(type);
}
// #ifdef SABINE
HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
{
HPREF_ELEMENT_TYPE type = HP_NONE;
int pnums[3];
int p[3];
INDEX_3 i3 (el.pnums[0], el.pnums[1], el.pnums[2]);
i3.Sort();
bool sing_face = faces.Used (i3);
// *testout << " facepoint " << facepoint << endl;
// Try all rotations of the trig
for (int j=0;j<3;j++)
{
int point_sing[3] = {0,0,0};
int edge_sing[3] = {0,0,0};
// *testout << " actual rotation of trig points " ;
for(int m=0;m<3;m++)
{
p[m] = (j+m)%3 +1; // local vertex number
pnums[m] = el.PNum(p[m]); // global vertex number
// *testout << pnums[m] << " \t ";
}
// *testout << endl ;
if(dim == 3)
{
// face point
for(int k=0;k<3;k++)
if(!sing_face)
{
// *testout << " fp [" << k << "] = " << facepoint[pnums[k]] << endl;
// *testout << " fd.DomainIn()" << fd.DomainIn() << endl;
// *testout << " fd.DomainOut()" << fd.DomainOut() << endl;
if( facepoint[pnums[k]] && (facepoint[pnums[k]] ==-1 ||
facepoint[pnums[k]] == fd.DomainIn() || facepoint[pnums[k]] == fd.DomainOut()))
point_sing[p[k]-1] = 1;
}
// if point is on face_edge in next step sing = 2
/* *testout << " pointsing NACH FACEPOints ... FALLS EDGEPOINT UMSETZEN" ;
for (int k=0;k<3;k++) *testout << "\t" << point_sing[p[k]-1] ;
*testout << endl; */
}
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1(TRIG);
if(dim==3)
{
for(int k=0;k<3;k++)
{
int ep1=p[eledges[k][0]-1];
int ep2=p[eledges[k][1]-1];
INDEX_2 i2(el.PNum(ep1),el.PNum(ep2));
if(edges.Used(i2))
{
edge_sing[k]=2;
point_sing[ep1-1] = 2;
point_sing[ep2-1] = 2;
}
else // face_edge?
{
i2.Sort();
if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1) // edge not face_edge acc. to surface in which trig lies
{
if(face_edges.Get(i2)==-1 ||face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut() )
{
edge_sing[k]=1;
}
else
{
point_sing[ep1-1] = 0; // set to edge_point
point_sing[ep2-1] = 0; // set to edge_point
}
}
}
/* *testout << " pointsing NACH edges UND FACEEDGES UMSETZEN ... " ;
for (int k=0;k<3;k++) *testout << "\t" << point_sing[p[k]-1] ;
*testout << endl;
*/
}
}
/*
*testout << " dim " << dim << endl;
*testout << " edgepoint_dom " << edgepoint_dom << endl;
*/
if(dim==2)
{
for(int k=0;k<3;k++)
{
int ep1=p[eledges[k][0]-1];
int ep2=p[eledges[k][1]-1];
INDEX_2 i2 = INDEX_2::Sort(el.PNum(ep1),el.PNum(ep2));
if(edges.Used(i2))
{
if(edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep1-1])) ||
edgepoint_dom.Used(INDEX_2(-1,pnums[ep1-1])) ||
edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep2-1])) ||
edgepoint_dom.Used(INDEX_2(-1,pnums[ep2-1])))
{
edge_sing[k]=2;
point_sing[ep1-1] = 2;
point_sing[ep2-1] = 2;
}
}
}
}
for (int k=0;k<3;k++)
if (edgepoint.Test(pnums[k]) &&
(dim==3 || edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[k])) || edgepoint_dom.Used(INDEX_2(-1,pnums[k]))))
//edgepoint, but not member of sing_edge on trig -> cp
{
INDEX_2 i2a=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+1)%3]));
INDEX_2 i2b=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+2)%3]));
if(!edges.Used(i2a) && !edges.Used(i2b))
point_sing[p[k]-1] = 3;
}
for(int k=0;k<3;k++)
if(cornerpoint.Test(el.PNum(p[k])))
point_sing[p[k]-1] = 3;
// *testout << "point_sing = " << point_sing[0] << point_sing[1] << point_sing[2] << endl;
if(edge_sing[0] + edge_sing[1] + edge_sing[2] == 0)
{
int ps = point_sing[0] + point_sing[1] + point_sing[2];
if(ps==0)
type = HP_TRIG;
else if(point_sing[p[0]-1] && !point_sing[p[1]-1] && !point_sing[p[2]-1])
type = HP_TRIG_SINGCORNER;
else if(point_sing[p[0]-1] && point_sing[p[1]-1] && !point_sing[p[2]-1])
type = HP_TRIG_SINGCORNER12;
else if(point_sing[p[0]-1] && point_sing[p[1]-1] && point_sing[p[2]-1])
{
if(dim==2) type = HP_TRIG_SINGCORNER123_2D;
else type = HP_TRIG_SINGCORNER123;
}
}
else
if (edge_sing[2] && !edge_sing[0] && !edge_sing[1]) //E[2]=(1,2)
{
int code = 0;
if(point_sing[p[0]-1] > edge_sing[2]) code+=1;
if(point_sing[p[1]-1] > edge_sing[2]) code+=2;
if(point_sing[p[2]-1]) code+=4;
HPREF_ELEMENT_TYPE types[] =
{
HP_TRIG_SINGEDGE,
HP_TRIG_SINGEDGECORNER1,
HP_TRIG_SINGEDGECORNER2,
HP_TRIG_SINGEDGECORNER12,
HP_TRIG_SINGEDGECORNER3,
HP_TRIG_SINGEDGECORNER13,
HP_TRIG_SINGEDGECORNER23,
HP_TRIG_SINGEDGECORNER123,
};
type = types[code];
} // E[0] = [0,2], E[1] =[1,2], E[2] = [0,1]
else
if(edge_sing[2] && !edge_sing[1] && edge_sing[0])
{
if(point_sing[p[2]-1] <= edge_sing[0] )
{
if(point_sing[p[1]-1]<= edge_sing[2]) type = HP_TRIG_SINGEDGES;
else type = HP_TRIG_SINGEDGES2;
}
else
{
if(point_sing[p[1]-1]<= edge_sing[2])
type = HP_TRIG_SINGEDGES3;
else type = HP_TRIG_SINGEDGES23;
}
}
else if (edge_sing[2] && edge_sing[1] && edge_sing[0])
type = HP_TRIG_3SINGEDGES;
// cout << " run for " << j << " gives type " << type << endl;
//*testout << " run for " << j << " gives type " << type << endl;
if(type!=HP_NONE) break;
}
// *testout << "type = " << type << endl;
for(int k=0;k<3;k++) el[k] = pnums[k];
/*if(type != HP_NONE)
{
cout << " TRIG with pnums " << pnums[0] << "\t" <<
pnums[1] << "\t" << pnums[2] << endl;
cout << " type " << type << endl;
}
*/
return(type);
}
#ifdef HPREF_OLD
HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
{
HPREF_ELEMENT_TYPE type = HP_NONE;
int pnums[3];
INDEX_3 i3 (el.pnums[0], el.pnums[1], el.pnums[2]);
i3.Sort();
bool sing_face = faces.Used (i3);
for (int j = 1; j <= 3; j++)
{
int ep1 = edgepoint.Test (el.PNumMod (j));
int ep2 = edgepoint.Test (el.PNumMod (j+1));
int ep3 = edgepoint.Test (el.PNumMod (j+2));
if (dim == 2)
{
// JS, Dec 11
ep1 = edgepoint_dom.Used (INDEX_2 (fd.SurfNr(), el.PNumMod(j))) ||
edgepoint_dom.Used (INDEX_2 (-1, el.PNumMod(j)));
ep2 = edgepoint_dom.Used (INDEX_2 (fd.SurfNr(), el.PNumMod(j+1))) ||
edgepoint_dom.Used (INDEX_2 (-1, el.PNumMod(j+1)));
ep3 = edgepoint_dom.Used (INDEX_2 (fd.SurfNr(), el.PNumMod(j+2))) ||
edgepoint_dom.Used (INDEX_2 (-1, el.PNumMod(j+2)));
/*
ep1 = edgepoint_dom.Used (INDEX_2 (el.index, el.PNumMod(j)));
ep2 = edgepoint_dom.Used (INDEX_2 (el.index, el.PNumMod(j+1)));
ep3 = edgepoint_dom.Used (INDEX_2 (el.index, el.PNumMod(j+2)));
*/
// ep3 = edgepoint_dom.Used (INDEX_2 (mesh.SurfaceElement(i).GetIndex(), el.PNumMod(j+2)));
}
int cp1 = cornerpoint.Test (el.PNumMod (j));
int cp2 = cornerpoint.Test (el.PNumMod (j+1));
int cp3 = cornerpoint.Test (el.PNumMod (j+2));
ep1 |= cp1;
ep2 |= cp2;
ep3 |= cp3;
// (*testout) << "cp = " << cp1 << cp2 << cp3 << ", ep = " << ep1 << ep2 << ep3 << endl;
int p[3] = { el.PNumMod (j), el.PNumMod (j+1), el.PNumMod (j+2)};
if(ep1)
{
INDEX_2 i2a=INDEX_2::Sort(p[0], p[1]);
INDEX_2 i2b=INDEX_2::Sort(p[0], p[2]);
if(!edges.Used(i2a) && !edges.Used(i2b))
cp1 = 1;
}
if(ep2)
{
INDEX_2 i2a=INDEX_2::Sort(p[1], p[0]);
INDEX_2 i2b=INDEX_2::Sort(p[1], p[2]);
if(!edges.Used(i2a) && !edges.Used(i2b))
cp2 = 1;
}
if(ep3)
{
INDEX_2 i2a=INDEX_2::Sort(p[2], p[0]);
INDEX_2 i2b=INDEX_2::Sort(p[2], p[1]);
if(!edges.Used(i2a) && !edges.Used(i2b))
cp3= 1;
}
int isedge1=0, isedge2=0, isedge3=0;
if(dim == 3 )
{
INDEX_2 i2;
i2 = INDEX_2(el.PNumMod (j), el.PNumMod (j+1));
isedge1 = edges.Used (i2);
i2.Sort();
if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1 &&
(face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) )
{
isedge1=1;
ep1 = 1; ep2=1;
}
i2 = INDEX_2(el.PNumMod (j+1), el.PNumMod (j+2));
isedge2 = edges.Used (i2);
i2.Sort();
if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1 &&
(face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) )
{
isedge2=1;
ep2 = 1; ep3=1;
}
i2 = INDEX_2(el.PNumMod (j+2), el.PNumMod (j+3));
isedge3 = edges.Used (i2);
i2.Sort();
if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1 &&
(face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) )
{
isedge3=1;
ep1 = 1; ep3=1;
}
// cout << " isedge " << isedge1 << " \t " << isedge2 << " \t " << isedge3 << endl;
if (!sing_face)
{
/*
if (!isedge1) { cp1 |= ep1; cp2 |= ep2; }
if (!isedge2) { cp2 |= ep2; cp3 |= ep3; }
if (!isedge3) { cp3 |= ep3; cp1 |= ep1; }
*/
ep1 |= facepoint [el.PNumMod(j)] != 0;
ep2 |= facepoint [el.PNumMod(j+1)] != 0;
ep3 |= facepoint [el.PNumMod(j+2)] != 0;
isedge1 |= face_edges.Used (INDEX_2::Sort (el.PNumMod(j), el.PNumMod(j+1)));
isedge2 |= face_edges.Used (INDEX_2::Sort (el.PNumMod(j+1), el.PNumMod(j+2)));
isedge3 |= face_edges.Used (INDEX_2::Sort (el.PNumMod(j+2), el.PNumMod(j+3)));
}
}
if(dim ==2)
{
INDEX_2 i2;
i2 = INDEX_2(el.PNumMod (j), el.PNumMod (j+1));
i2.Sort();
isedge1 = edges.Used (i2);
if(isedge1)
{
ep1 = 1; ep2=1;
}
i2 = INDEX_2(el.PNumMod (j+1), el.PNumMod (j+2));
i2.Sort();
isedge2 = edges.Used (i2);
if(isedge2)
{
ep2 = 1; ep3=1;
}
i2 = INDEX_2(el.PNumMod (j+2), el.PNumMod (j+3));
i2.Sort();
isedge3 = edges.Used (i2);
if(isedge3)
{
ep1 = 1; ep3=1;
}
}
/*
cout << " used " << face_edges.Used (INDEX_2::Sort (el.PNumMod(j), el.PNumMod(j+1))) << endl;
cout << " isedge " << isedge1 << " \t " << isedge2 << " \t " << isedge3 << endl;
cout << " ep " << ep1 << "\t" << ep2 << " \t " << ep3 << endl;
cout << " cp " << cp1 << "\t" << cp2 << " \t " << cp3 << endl;
*/
if (isedge1 + isedge2 + isedge3 == 0)
{
if (!ep1 && !ep2 && !ep3)
type = HP_TRIG;
if (ep1 && !ep2 && !ep3)
type = HP_TRIG_SINGCORNER;
if (ep1 && ep2 && !ep3)
type = HP_TRIG_SINGCORNER12;
if (ep1 && ep2 && ep3)
{
if (dim == 2)
type = HP_TRIG_SINGCORNER123_2D;
else
type = HP_TRIG_SINGCORNER123;
}
if (type != HP_NONE)
{
pnums[0] = el.PNumMod (j);
pnums[1] = el.PNumMod (j+1);
pnums[2] = el.PNumMod (j+2);
break;
}
}
if (isedge1 && !isedge2 && !isedge3)
{
int code = 0;
if (cp1) code += 1;
if (cp2) code += 2;
if (ep3) code += 4;
HPREF_ELEMENT_TYPE types[] =
{
HP_TRIG_SINGEDGE,
HP_TRIG_SINGEDGECORNER1,
HP_TRIG_SINGEDGECORNER2,
HP_TRIG_SINGEDGECORNER12,
HP_TRIG_SINGEDGECORNER3,
HP_TRIG_SINGEDGECORNER13,
HP_TRIG_SINGEDGECORNER23,
HP_TRIG_SINGEDGECORNER123,
};
type = types[code];
pnums[0] = el.PNumMod (j);
pnums[1] = el.PNumMod (j+1);
pnums[2] = el.PNumMod (j+2);
break;
}
if (isedge1 && !isedge2 && isedge3)
{
if (!cp3)
{
if (!cp2) type = HP_TRIG_SINGEDGES;
else type = HP_TRIG_SINGEDGES2;
}
else
{
if (!cp2) type = HP_TRIG_SINGEDGES3;
else type = HP_TRIG_SINGEDGES23;
}
pnums[0] = el.PNumMod (j);
pnums[1] = el.PNumMod (j+1);
pnums[2] = el.PNumMod (j+2);
break;
}
if (isedge1 && isedge2 && isedge3)
{
type = HP_TRIG_3SINGEDGES;
pnums[0] = el.PNumMod (j);
pnums[1] = el.PNumMod (j+1);
pnums[2] = el.PNumMod (j+2);
break;
}
}
for(int k=0;k<3;k++) el[k] = pnums[k];
/*if(type != HP_NONE)
{
cout << " TRIG with pnums " << pnums[0] << "\t" <<
pnums[1] << "\t" << pnums[2] << endl;
cout << " type " << type << endl;
}
*/
return(type);
}
#endif
HPREF_ELEMENT_TYPE ClassifyQuad(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
{
HPREF_ELEMENT_TYPE type = HP_NONE;
int ep1(-1), ep2(-1), ep3(-1), ep4(-1), cp1(-1), cp2(-1), cp3(-1), cp4(-1);
int isedge1, isedge2, isedge3, isedge4;
// *testout << "edges = " << edges << endl;
for (int j = 1; j <= 4; j++)
{
ep1 = edgepoint.Test (el.PNumMod (j));
ep2 = edgepoint.Test (el.PNumMod (j+1));
ep3 = edgepoint.Test (el.PNumMod (j+2));
ep4 = edgepoint.Test (el.PNumMod (j+3));
if (dim == 2)
{
ep1 = edgepoint_dom.Used (INDEX_2 (el.GetIndex(), el.PNumMod(j)));
ep2 = edgepoint_dom.Used (INDEX_2 (el.GetIndex(), el.PNumMod(j+1)));
ep3 = edgepoint_dom.Used (INDEX_2 (el.GetIndex(), el.PNumMod(j+2)));
ep4 = edgepoint_dom.Used (INDEX_2 (el.GetIndex(), el.PNumMod(j+3)));
}
cp1 = cornerpoint.Test (el.PNumMod (j));
cp2 = cornerpoint.Test (el.PNumMod (j+1));
cp3 = cornerpoint.Test (el.PNumMod (j+2));
cp4 = cornerpoint.Test (el.PNumMod (j+3));
ep1 |= cp1;
ep2 |= cp2;
ep3 |= cp3;
ep4 |= cp4;
int p[4] = { el.PNumMod (j), el.PNumMod (j+1), el.PNumMod (j+2), el.PNumMod(j+4)};
//int epp[4] = { ep1, ep2, ep3, ep4};
int cpp[4] = { cp1, cp2, cp3, cp4};
for(int k=0;k<0;k++)
{
INDEX_2 i2a=INDEX_2::Sort(p[k], p[(k+1)%4]);
INDEX_2 i2b=INDEX_2::Sort(p[k], p[(k-1)%4]);
if(!edges.Used(i2a) && !edges.Used(i2b))
cpp[k] = 1;
}
cp1= cpp[0]; cp2=cpp[1]; cp3=cpp[2]; cp4=cpp[3];
if(dim ==3)
{
INDEX_2 i2;
i2 = INDEX_2(el.PNumMod (j), el.PNumMod (j+1));
// i2.Sort();
isedge1 = edges.Used (i2);
i2.Sort();
if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1 &&
(face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) )
{
isedge1=1;
ep1 = 1; ep2=1;
}
i2 = INDEX_2(el.PNumMod (j+1), el.PNumMod (j+2));
// i2.Sort();
isedge2 = edges.Used (i2);
i2.Sort();
if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1 &&
(face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) )
{
isedge2=1;
ep2=1; ep3=1;
}
i2 = INDEX_2(el.PNumMod (j+2), el.PNumMod (j+3));
// i2.Sort();
isedge3 = edges.Used (i2);
i2.Sort();
if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1 && (face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) )
{
isedge3=1;
ep3=1; ep4=1;
}
i2 = INDEX_2(el.PNumMod (j+3), el.PNumMod (j+4));
// i2.Sort();
isedge4 = edges.Used (i2);
i2.Sort();
if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1 &&
(face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) )
{
isedge4=1;
ep4=1; ep1=1;
}
//MH***********************************************************************************************************
if(ep1)
if(edgepoint.Test(p[0]))
{
INDEX_2 i2a=INDEX_2::Sort(p[0], p[1]);
INDEX_2 i2b=INDEX_2::Sort(p[0], p[3]);
if(!edges.Used(i2a) && !edges.Used(i2b))
cp1 = 1;
}
if(ep2)
if(edgepoint.Test(p[1]))
{
INDEX_2 i2a=INDEX_2::Sort(p[0], p[1]);
INDEX_2 i2b=INDEX_2::Sort(p[1], p[2]);
if(!edges.Used(i2a) && !edges.Used(i2b))
cp2 = 1;
}
if(ep3)
if(edgepoint.Test(p[2]))
{
INDEX_2 i2a=INDEX_2::Sort(p[2], p[1]);
INDEX_2 i2b=INDEX_2::Sort(p[3], p[2]);
if(!edges.Used(i2a) && !edges.Used(i2b))
cp3 = 1;
}
if(ep4)
if(edgepoint.Test(p[3]))
{
INDEX_2 i2a=INDEX_2::Sort(p[0], p[3]);
INDEX_2 i2b=INDEX_2::Sort(p[3], p[2]);
if(!edges.Used(i2a) && !edges.Used(i2b))
cp4 = 1;
}
//MH*****************************************************************************************************************************
}
else
{
INDEX_2 i2;
i2 = INDEX_2(el.PNumMod (j), el.PNumMod (j+1));
i2.Sort();
isedge1 = edges.Used (i2);
if(isedge1)
{
ep1 = 1; ep2=1;
}
i2 = INDEX_2(el.PNumMod (j+1), el.PNumMod (j+2));
i2.Sort();
isedge2 = edges.Used (i2);
if(isedge2)
{
ep2=1; ep3=1;
}
i2 = INDEX_2(el.PNumMod (j+2), el.PNumMod (j+3));
i2.Sort();
isedge3 = edges.Used (i2);
if(isedge3)
{
ep3=1; ep4=1;
}
i2 = INDEX_2(el.PNumMod (j+3), el.PNumMod (j+4));
i2.Sort();
isedge4 = edges.Used (i2);
if(isedge4)
{
ep4=1; ep1=1;
}
}
int sumcp = cp1 + cp2 + cp3 + cp4;
int sumep = ep1 + ep2 + ep3 + ep4;
int sumedge = isedge1 + isedge2 + isedge3 + isedge4;
*testout << "isedge = " << isedge1 << isedge2 << isedge3 << isedge4 << endl;
*testout << "iscp = " << cp1 << cp2 << cp3 << cp4 << endl;
*testout << "isep = " << ep1 << ep2 << ep3 << ep4 << endl;
switch (sumedge)
{
case 0:
{
switch (sumep)
{
case 0:
type = HP_QUAD;
break;
case 1:
if (ep1) type = HP_QUAD_SINGCORNER;
break;
case 2:
{
if (ep1 && ep2) type = HP_QUAD_0E_2VA;
if (ep1 && ep3) type = HP_QUAD_0E_2VB;
break;
}
case 3:
if (!ep4) type = HP_QUAD_0E_3V;
break;
case 4:
type = HP_QUAD_0E_4V;
break;
}
break;
}
case 1:
{
if (isedge1)
{
switch (cp1+cp2+ep3+ep4)
{
case 0:
type = HP_QUAD_SINGEDGE;
break;
case 1:
{
if (cp1) type = HP_QUAD_1E_1VA;
if (cp2) type = HP_QUAD_1E_1VB;
if (ep3) type = HP_QUAD_1E_1VC;
if (ep4) type = HP_QUAD_1E_1VD;
break;
}
case 2:
{
if (cp1 && cp2) type = HP_QUAD_1E_2VA;
if (cp1 && ep3) type = HP_QUAD_1E_2VB;
if (cp1 && ep4) type = HP_QUAD_1E_2VC;
if (cp2 && ep3) type = HP_QUAD_1E_2VD;
if (cp2 && ep4) type = HP_QUAD_1E_2VE;
if (ep3 && ep4) type = HP_QUAD_1E_2VF;
break;
}
case 3:
{
if (cp1 && cp2 && ep3) type = HP_QUAD_1E_3VA;
if (cp1 && cp2 && ep4) type = HP_QUAD_1E_3VB;
if (cp1 && ep3 && ep4) type = HP_QUAD_1E_3VC;
if (cp2 && ep3 && ep4) type = HP_QUAD_1E_3VD;
break;
}
case 4:
{
type = HP_QUAD_1E_4V;
break;
}
}
}
break;
}
case 2:
{
if (isedge1 && isedge4)
{
if (!cp2 && !ep3 && !cp4)
type = HP_QUAD_2E;
if (cp2 && !ep3 && !cp4)
type = HP_QUAD_2E_1VA;
if (!cp2 && ep3 && !cp4)
type = HP_QUAD_2E_1VB;
if (!cp2 && !ep3 && cp4)
type = HP_QUAD_2E_1VC;
if (cp2 && ep3 && !cp4)
type = HP_QUAD_2E_2VA;
if (cp2 && !ep3 && cp4)
type = HP_QUAD_2E_2VB;
if (!cp2 && ep3 && cp4)
type = HP_QUAD_2E_2VC;
if (cp2 && ep3 && cp4)
type = HP_QUAD_2E_3V;
}
if (isedge1 && isedge3)
{
switch (sumcp)
{
case 0:
type = HP_QUAD_2EB_0V; break;
case 1:
{
if (cp1) type = HP_QUAD_2EB_1VA;
if (cp2) type = HP_QUAD_2EB_1VB;
break;
}
case 2:
{
if (cp1 && cp2) { type = HP_QUAD_2EB_2VA; }
if (cp1 && cp3) { type = HP_QUAD_2EB_2VB; }
if (cp1 && cp4) { type = HP_QUAD_2EB_2VC; }
if (cp2 && cp4) { type = HP_QUAD_2EB_2VD; }
break;
}
case 3:
{
if (cp1 && cp2 && cp3) { type = HP_QUAD_2EB_3VA; }
if (cp1 && cp2 && cp4) { type = HP_QUAD_2EB_3VB; }
break;
}
case 4:
{
type = HP_QUAD_2EB_4V; break;
}
}
}
break;
}
case 3:
{
if (isedge1 && isedge2 && isedge4)
{
if (!cp3 && !cp4) type = HP_QUAD_3E;
if (cp3 && !cp4) type = HP_QUAD_3E_3VA;
if (!cp3 && cp4) type = HP_QUAD_3E_3VB;
if (cp3 && cp4) type = HP_QUAD_3E_4V;
}
break;
}
case 4:
{
type = HP_QUAD_4E;
break;
}
}
if (type != HP_NONE)
{
int pnums[4];
pnums[0] = el.PNumMod (j);
pnums[1] = el.PNumMod (j+1);
pnums[2] = el.PNumMod (j+2);
pnums[3] = el.PNumMod (j+3);
for (int k=0;k<4;k++) el[k] = pnums[k];
/* cout << " QUAD with pnums " << pnums[0] << "\t" <<
pnums[1] << "\t" << pnums[2] << "\t" << pnums[3]
<< endl << " of type " << type << endl; */
break;
}
}
if (type == HP_NONE)
{
(*testout) << "undefined element" << endl
<< "cp = " << cp1 << cp2 << cp3 << cp4 << endl
<< "ep = " << ep1 << ep2 << ep3 << ep4 << endl
<< "isedge = " << isedge1 << isedge2 << isedge3
<< isedge4 << endl;
}
*testout << "quad type = " << type << endl;
return type;
}
HPREF_ELEMENT_TYPE ClassifyHex(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
{
HPREF_ELEMENT_TYPE type = HP_NONE;
// implementation only for HP_HEX_1F_0E_0V
// HP_HEX_1FA_1FB_0E_0V
// HP_HEX
// up to now other cases are refine dummies
// indices of bot,top-faces combinations
int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}};
int p[8];
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (HEX);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (HEX);
for(int m=0;m<6 && type == HP_NONE;m++)
for(int j=0;j<4 && type == HP_NONE;j++)
{
int point_sing[8]={0,0,0,0,0,0,0,0};
int face_sing[6] = {0,0,0,0,0,0};
int edge_sing[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
int spoint=0, sface=0, sedge=0;
for(int l=0;l<4;l++)
{
p[l] = elfaces[index[m][0]][(4-j-l)%4];
p[l+4] = elfaces[index[m][1]][(j+l)%4];
}
for(int l=0;l<8;l++)
if(cornerpoint.Test(el.PNum(p[l])))
{
point_sing[p[l]-1]=3;
spoint++;
}
else if(edgepoint.Test(el.PNum(p[l]))) point_sing[p[l]-1]=2;
else if (facepoint[el.PNum(p[l])] == -1 || facepoint[el.PNum(p[l])] == el.GetIndex())
point_sing[p[l]-1] = 1;
for(int k=0;k<12;k++)
{
INDEX_2 i2 = INDEX_2 :: Sort(el.PNum(p[eledges[k][0]-1]),el.PNum(p[eledges[k][1]-1]));
if (edges.Used(i2))
{
edge_sing[k] = 2;
sedge++;
}
else edge_sing[k] = face_edges.Used(i2);
}
for (int k=0;k<6;k++)
{
INDEX_3 i3;
INDEX_4 i4 = INDEX_4(el.pnums[p[elfaces[k][0]-1]-1], el.pnums[p[elfaces[k][1]-1]-1], el.pnums[p[elfaces[k][2]-1]-1],el.pnums[p[elfaces[k][3]-1]-1]);
i4.Sort();
i3 = INDEX_3(i4.I1(), i4.I2(), i4.I3());
if (faces.Used (i3))
{
int domnr = faces.Get(i3);
if (domnr == -1 || domnr == el.GetIndex())
{
face_sing[k] = 1;
sface++;
}
}
}
if(!sface && !sedge && !spoint) type = HP_HEX;
if(!sedge && !spoint)
{
if(face_sing[0] && face_sing[2] && sface==2)
type = HP_HEX_1FA_1FB_0E_0V;
if (face_sing[0] && sface==1)
type = HP_HEX_1F_0E_0V;
}
el.type=type;
if(type != HP_NONE)
{
int pnums[8];
for(int l=0;l<8;l++) pnums[l] = el[p[l]-1];
for(int l=0;l<8;l++) el[l] = pnums[l];
/* cout << " HEX with pnums " << pnums[0] << "\t" <<
pnums[1] << "\t" << pnums[2] << "\t" << pnums[3] << "\t" <<
pnums[4] << "\t" << pnums[5] << endl << " of type " << type << endl; */
break;
}
}
return (type);
}
HPREF_ELEMENT_TYPE ClassifyHex7 (HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
{
// HPREF_ELEMENT_TYPE type = HP_NONE;
// no singular
// singular bottom
// singular top
// indices of bot,top-faces combinations
// int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}};
// int p[8];
// const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (HEX);
// const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (HEX);
INDEX_4 fbot4 = { el.pnums[0], el.pnums[1], el.pnums[2], el.pnums[3] };
INDEX_3 ftop = { el.pnums[4], el.pnums[5], el.pnums[6] };
fbot4.Sort();
INDEX_3 fbot = { fbot4[0], fbot4[1], fbot4[2] };
ftop.Sort();
bool singbot = faces.Used(fbot);
bool singtop = faces.Used(ftop);
if (singbot)
el.type = HP_HEX7_1FA;
else if (singtop)
el.type = HP_HEX7_1FB;
else
el.type = HP_HEX7;
return el.type;
}
HPREF_ELEMENT_TYPE ClassifySegm(HPRefElement & hpel, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
{
int cp1 = cornerpoint.Test (hpel[0]);
int cp2 = cornerpoint.Test (hpel[1]);
INDEX_2 i2;
i2 = INDEX_2(hpel[0], hpel[1]);
i2.Sort();
if (!edges.Used (i2))
{
cp1 = edgepoint.Test (hpel[0]);
cp2 = edgepoint.Test (hpel[1]);
}
if(!edges.Used(i2) && !face_edges.Used(i2))
{
if(facepoint[hpel[0]]!=0) cp1=1;
if(facepoint[hpel[1]]!=0) cp2=1;
}
if(edges.Used(i2) && !face_edges.Used(i2))
{
if(facepoint[hpel[0]]) cp1 = 1;
if(facepoint[hpel[1]]) cp2 = 1;
}
if (!cp1 && !cp2)
hpel.type = HP_SEGM;
else if (cp1 && !cp2)
hpel.type = HP_SEGM_SINGCORNERL;
else if (!cp1 && cp2)
hpel.type = HP_SEGM_SINGCORNERR;
else
hpel.type = HP_SEGM_SINGCORNERS;
// cout << " SEGM found with " << hpel[0] << " \t" << hpel[1] << endl << " of type " << hpel.type << endl;
return(hpel.type) ;
}
HPREF_ELEMENT_TYPE ClassifyPyramid(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
{
// *testout << "classify pyramid, pnums = ";
// for (int i = 0; i < 5; i++) *testout << el.pnums[i] << " ";
// *testout << endl;
HPREF_ELEMENT_TYPE type = HP_NONE;
// implementation only for HP_PYRAMID
// HP_PYRAMID_0E_1V
// HP_PYRAMID_EDGES
// HP_PYRAMID_1FB_0E_1VA
// up to now other cases are refine dummies
// indices of bot,top-faces combinations
// int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}};
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (PYRAMID);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (PYRAMID);
int point_sing[5]={0,0,0,0,0};
int face_sing[5] = {0,0,0,0,0};
int edge_sing[8] = {0,0,0,0,0,0,0,0};
int spoint=0, sedge=0, sface=0;
for(int m=0;m<4 && type == HP_NONE;m++)
{
*testout << "m = " << m << endl;
int p[5] = {m%4, m%4+1, m%4+2, m%4+3, 4};
for(int l=0;l<5;l++)
{
if(cornerpoint.Test(el.pnums[p[l]]))
point_sing[l]=3;
else if(edgepoint.Test(el.pnums[p[l]]))
point_sing[l]=2;
else if (facepoint[el.pnums[p[l]]] == -1 || facepoint[el.pnums[p[l]]] == el.GetIndex())
point_sing[l] = 1;
spoint += point_sing[l];
}
for(int k=0;k<8;k++)
{
INDEX_2 i2 = INDEX_2 :: Sort(el.pnums[p[eledges[k][0]-1]],
el.pnums[p[eledges[k][1]-1]]);
if (edges.Used(i2))
edge_sing[k] = 2;
else
edge_sing[k] = face_edges.Used(i2);
sedge += edge_sing[k];
}
for (int k=0;k<5;k++)
{
INDEX_3 i3;
/*
INDEX_4 i4 = INDEX_4(el.pnums[p[elfaces[k][0]-1]], el.pnums[p[elfaces[k][1]-1]], el.pnums[p[elfaces[k][2]-1]],
el.pnums[p[elfaces[k][3]-1]]);
i4.Sort();
i3 = INDEX_3(i4.I1(), i4.I2(), i4.I3());
*/
if (k < 4)
{
i3 = INDEX_3(el.pnums[p[elfaces[k][0]-1]], el.pnums[p[elfaces[k][1]-1]], el.pnums[p[elfaces[k][2]-1]]);
i3.Sort();
}
else
{
INDEX_4 i4 = INDEX_4(el.pnums[p[elfaces[k][0]-1]], el.pnums[p[elfaces[k][1]-1]], el.pnums[p[elfaces[k][2]-1]],
el.pnums[p[elfaces[k][3]-1]]);
i4.Sort();
i3 = INDEX_3(i4.I1(), i4.I2(), i4.I3());
}
if (faces.Used (i3))
{
int domnr = faces.Get(i3);
if (domnr == -1 || domnr == el.GetIndex())
face_sing[k] = 1;
}
sface +=face_sing[k];
}
*testout << "point_sing: ";
for (int k = 0; k < 5; k++) *testout << point_sing[k] << " ";
*testout << endl;
*testout << "edge_sing: ";
for (int k = 0; k < 4; k++) *testout << edge_sing[k] << " ";
*testout << endl;
*testout << "face_sing: ";
for (int k = 0; k < 5; k++) *testout << face_sing[k] << " ";
*testout << endl;
if(!sface && !spoint && !sedge) return(HP_PYRAMID);
if(!sface && !sedge && point_sing[p[0]] == spoint)
type = HP_PYRAMID_0E_1V;
if(!sface && edge_sing[0] + edge_sing[2] == sedge &&
spoint == point_sing[0] + point_sing[1] + point_sing[3])
type = HP_PYRAMID_EDGES;
if(sface && sface == face_sing[0] && spoint == point_sing[4] + 2)
{
if (point_sing[4] == 1)
type = HP_PYRAMID_1FB_0E_0V;
else
type = HP_PYRAMID_1FB_0E_1VA;
}
if(type != HP_NONE)
{
int pnums[8];
for(int l=0;l<5;l++) pnums[l] = el[p[l]];
for(int l=0;l<5;l++) el[l] = pnums[l];
el.type=type;
break;
}
}
return (type);
}