ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-04-12 16:17:29 +05:00
Code Issues Packages Projects Releases Wiki Activity

reduced memory requirement for the setup of topology faces

Browse Source
This commit is contained in:
Joachim Schoeberl 2009-10-27 18:14:45 +00:00
parent aa805dfafe
commit 6f39164242
10 changed files with 2054 additions and 1473 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
libsrc/general/hashtabl.hpp
View File

@ -717,6 +717,27 @@ public:
} }
}; };
template <typename T>
inline ostream & operator<< (ostream & ost, const INDEX_2_CLOSED_HASHTABLE<T> & ht)
{
for (int i = 0; i < ht.Size(); i++)
if (ht.UsedPos(i))
{
INDEX_2 hash;
T data;
ht.GetData (i, hash, data);
ost << "hash = " << hash << ", data = " << data << endl;
}
return ost;
}
class BASE_INDEX_3_CLOSED_HASHTABLE class BASE_INDEX_3_CLOSED_HASHTABLE
{ {
protected: protected:
@ -724,7 +745,6 @@ protected:
int invalid; int invalid;
protected: protected:
// public: //SZ
BASE_INDEX_3_CLOSED_HASHTABLE (int size) BASE_INDEX_3_CLOSED_HASHTABLE (int size)
: hash(size) : hash(size)
{ {
@ -896,6 +916,20 @@ public:
template <typename T>
inline ostream & operator<< (ostream & ost, const INDEX_3_CLOSED_HASHTABLE<T> & ht)
{
for (int i = 0; i < ht.Size(); i++)
if (ht.UsedPos(i))
{
INDEX_3 hash;
T data;
ht.GetData (i, hash, data);
ost << "hash = " << hash << ", data = " << data << endl;
}
return ost;
}

7
libsrc/general/profiler.cpp
View File

@ -42,8 +42,8 @@ namespace netgen
// thus we use the C-variant: // thus we use the C-variant:
if (getenv ("NGPROFILE"))
/* {
char filename[100]; char filename[100];
#ifdef PARALLEL #ifdef PARALLEL
sprintf (filename, "netgen.prof.%d", id); sprintf (filename, "netgen.prof.%d", id);
@ -51,10 +51,11 @@ namespace netgen
sprintf (filename, "netgen.prof"); sprintf (filename, "netgen.prof");
#endif #endif
printf ("write profile to file %s\n", filename);
FILE *prof = fopen(filename,"w"); FILE *prof = fopen(filename,"w");
Print (prof); Print (prof);
fclose(prof); fclose(prof);
*/ }
} }

18
libsrc/meshing/bisect.cpp
View File

@ -2643,6 +2643,10 @@ namespace netgen
PrintMessage(1,"Mesh bisection"); PrintMessage(1,"Mesh bisection");
PushStatus("Mesh bisection"); PushStatus("Mesh bisection");
static int timer = NgProfiler::CreateTimer ("Bisect");
NgProfiler::RegionTimer reg1 (timer);
static int localizetimer = NgProfiler::CreateTimer("localize edgepoints"); static int localizetimer = NgProfiler::CreateTimer("localize edgepoints");
NgProfiler::RegionTimer * loct = new NgProfiler::RegionTimer(localizetimer); NgProfiler::RegionTimer * loct = new NgProfiler::RegionTimer(localizetimer);
@ -3524,6 +3528,7 @@ namespace netgen
} }
while (hangingvol || hangingsurf || hangingedge); while (hangingvol || hangingsurf || hangingedge);
/*
if (printmessage_importance>0) if (printmessage_importance>0)
{ {
ostringstream strstr; ostringstream strstr;
@ -3536,8 +3541,16 @@ namespace netgen
} }
strstr << mesh.GetNP() << " points"; strstr << mesh.GetNP() << " points";
PrintMessage(4,strstr.str()); PrintMessage(4,strstr.str());
} }
*/
PrintMessage (4, mtets.Size(), " tets");
PrintMessage (4, mtris.Size(), " trigs");
if (mprisms.Size())
{
PrintMessage (4, mprisms.Size(), " prisms");
PrintMessage (4, mquads.Size(), " quads");
}
PrintMessage (4, mesh.GetNP(), " points");
} }
@ -3859,9 +3872,6 @@ namespace netgen
// Check/Repair // Check/Repair
//cout << "Hallo Welt" << endl;
//getchar();
static bool repaired_once; static bool repaired_once;
if(mesh.mglevels == 1) if(mesh.mglevels == 1)
repaired_once = false; repaired_once = false;

14
libsrc/meshing/classifyhpel.hpp
View File

@ -455,7 +455,7 @@ HPREF_ELEMENT_TYPE ClassifyPrism(HPRefElement & el, INDEX_2_HASHTABLE<int> & edg
point_sing[p[j]-1] = 1; point_sing[p[j]-1] = 1;
} }
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (PRISM); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (PRISM);
for(int k=0;k<9;k++) 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])); INDEX_2 i2 = INDEX_2 :: Sort(el.PNum(p[eledges[k][0]-1]),el.PNum(p[eledges[k][1]-1]));
@ -463,7 +463,7 @@ HPREF_ELEMENT_TYPE ClassifyPrism(HPRefElement & el, INDEX_2_HASHTABLE<int> & edg
else edge_sing[k] = face_edges.Used(i2); else edge_sing[k] = face_edges.Used(i2);
} }
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces (PRISM); const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (PRISM);
for (int k=0;k<5;k++) for (int k=0;k<5;k++)
{ {
INDEX_3 i3; INDEX_3 i3;
@ -710,7 +710,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
*testout << endl; */ *testout << endl; */
} }
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges(TRIG); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1(TRIG);
if(dim==3) if(dim==3)
{ {
@ -1499,8 +1499,8 @@ HPREF_ELEMENT_TYPE ClassifyHex(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges
// indices of bot,top-faces combinations // indices of bot,top-faces combinations
int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}}; int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}};
int p[8]; int p[8];
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces (HEX); const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (HEX);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (HEX); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (HEX);
for(int m=0;m<6 && type == HP_NONE;m++) for(int m=0;m<6 && type == HP_NONE;m++)
for(int j=0;j<4 && type == HP_NONE;j++) for(int j=0;j<4 && type == HP_NONE;j++)
@ -1643,8 +1643,8 @@ HPREF_ELEMENT_TYPE ClassifyPyramid(HPRefElement & el, INDEX_2_HASHTABLE<int> & e
// indices of bot,top-faces combinations // indices of bot,top-faces combinations
// int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}}; // int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}};
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces (PYRAMID); const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (PYRAMID);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (PYRAMID); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (PYRAMID);
int point_sing[5]={0,0,0,0,0}; int point_sing[5]={0,0,0,0,0};
int face_sing[5] = {0,0,0,0,0}; int face_sing[5] = {0,0,0,0,0};

36
libsrc/meshing/curvedelems.cpp
View File

@ -545,7 +545,7 @@ namespace netgen
top.GetSurfaceElementEdges (i+1, edgenrs); top.GetSurfaceElementEdges (i+1, edgenrs);
for (int j = 0; j < edgenrs.Size(); j++) for (int j = 0; j < edgenrs.Size(); j++)
edgenrs[j]--; edgenrs[j]--;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (el.GetType()); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (el.GetType());
for (int i2 = 0; i2 < edgenrs.Size(); i2++) for (int i2 = 0; i2 < edgenrs.Size(); i2++)
{ {
@ -1371,7 +1371,7 @@ namespace netgen
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
shapes(j) = lami[j] * lami[j]; shapes(j) = lami[j] * lami[j];
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TRIG); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TRIG);
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
{ {
double wi = edgeweight[info.edgenrs[j]]; double wi = edgeweight[info.edgenrs[j]];
@ -1394,7 +1394,7 @@ namespace netgen
if (info.order == 1) return; if (info.order == 1) return;
int ii = 3; int ii = 3;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TRIG); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TRIG);
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
{ {
@ -1465,7 +1465,7 @@ namespace netgen
}; };
int ii = 4; int ii = 4;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (QUAD); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (QUAD);
for (int i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
{ {
@ -1525,7 +1525,7 @@ namespace netgen
dshapes(j,1) = 2 * lami[j] * dlami[j][1]; dshapes(j,1) = 2 * lami[j] * dlami[j][1];
} }
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TRIG); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TRIG);
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
{ {
double wi = edgeweight[info.edgenrs[j]]; double wi = edgeweight[info.edgenrs[j]];
@ -1573,7 +1573,7 @@ namespace netgen
lami[2] = 1-xi(0)-xi(1); lami[2] = 1-xi(0)-xi(1);
int ii = 3; int ii = 3;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TRIG); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TRIG);
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
{ {
@ -1707,7 +1707,7 @@ namespace netgen
ArrayMem<double, 20> hshapes(order+1), hdshapes(order+1); ArrayMem<double, 20> hshapes(order+1), hdshapes(order+1);
int ii = 4; int ii = 4;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (QUAD); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (QUAD);
for (int i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
{ {
@ -1987,7 +1987,7 @@ namespace netgen
for (int j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
shapes(j) = lami[j] * lami[j]; shapes(j) = lami[j] * lami[j];
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TET); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TET);
for (int j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
{ {
double wi = edgeweight[info.edgenrs[j]]; double wi = edgeweight[info.edgenrs[j]];
@ -2013,7 +2013,7 @@ namespace netgen
if (info.order == 1) return; if (info.order == 1) return;
int ii = 4; int ii = 4;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TET); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TET);
for (int i = 0; i < 6; i++) for (int i = 0; i < 6; i++)
{ {
int eorder = edgeorder[info.edgenrs[i]]; int eorder = edgeorder[info.edgenrs[i]];
@ -2026,7 +2026,7 @@ namespace netgen
ii += eorder-1; ii += eorder-1;
} }
} }
const ELEMENT_FACE * faces = MeshTopology::GetFaces (TET); const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (TET);
for (int i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
{ {
int forder = faceorder[info.facenrs[i]]; int forder = faceorder[info.facenrs[i]];
@ -2091,7 +2091,7 @@ namespace netgen
int ii = 6; int ii = 6;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (PRISM); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (PRISM);
for (int i = 0; i < 6; i++) // horizontal edges for (int i = 0; i < 6; i++) // horizontal edges
{ {
int eorder = edgeorder[info.edgenrs[i]]; int eorder = edgeorder[info.edgenrs[i]];
@ -2131,7 +2131,7 @@ namespace netgen
} }
// FACE SHAPES // FACE SHAPES
const ELEMENT_FACE * faces = MeshTopology::GetFaces (PRISM); const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (PRISM);
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
{ {
int forder = faceorder[info.facenrs[i]]; int forder = faceorder[info.facenrs[i]];
@ -2170,7 +2170,7 @@ namespace netgen
if (info.order == 1) return; if (info.order == 1) return;
int ii = 5; int ii = 5;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (PYRAMID); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (PYRAMID);
for (int i = 0; i < 4; i++) // horizontal edges for (int i = 0; i < 4; i++) // horizontal edges
{ {
int eorder = edgeorder[info.edgenrs[i]]; int eorder = edgeorder[info.edgenrs[i]];
@ -2245,7 +2245,7 @@ namespace netgen
dshapes(j,2) = 2 * lami[j] * dlami[j][2]; dshapes(j,2) = 2 * lami[j] * dlami[j][2];
} }
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TET); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TET);
for (int j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
{ {
double wi = edgeweight[info.edgenrs[j]]; double wi = edgeweight[info.edgenrs[j]];
@ -2283,7 +2283,7 @@ namespace netgen
double lami[] = { xi(0), xi(1), xi(2), 1-xi(0)-xi(1)-xi(2) }; double lami[] = { xi(0), xi(1), xi(2), 1-xi(0)-xi(1)-xi(2) };
int ii = 4; int ii = 4;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TET); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TET);
for (int i = 0; i < 6; i++) for (int i = 0; i < 6; i++)
{ {
int eorder = edgeorder[info.edgenrs[i]]; int eorder = edgeorder[info.edgenrs[i]];
@ -2314,7 +2314,7 @@ namespace netgen
} }
} }
const ELEMENT_FACE * faces = MeshTopology::GetFaces (TET); const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (TET);
for (int i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
{ {
int forder = faceorder[info.facenrs[i]]; int forder = faceorder[info.facenrs[i]];
@ -2428,7 +2428,7 @@ namespace netgen
if (info.order == 1) return; if (info.order == 1) return;
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (PRISM); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (PRISM);
for (int i = 0; i < 6; i++) // horizontal edges for (int i = 0; i < 6; i++) // horizontal edges
{ {
int order = edgeorder[info.edgenrs[i]]; int order = edgeorder[info.edgenrs[i]];
@ -2502,7 +2502,7 @@ namespace netgen
if (info.order == 2) return; if (info.order == 2) return;
// FACE SHAPES // FACE SHAPES
const ELEMENT_FACE * faces = MeshTopology::GetFaces (PRISM); const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (PRISM);
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
{ {
int forder = faceorder[info.facenrs[i]]; int forder = faceorder[info.facenrs[i]];

8
libsrc/meshing/hprefinement.cpp
View File

@ -991,7 +991,7 @@ namespace netgen
int npi = mesh.AddPoint (center); int npi = mesh.AddPoint (center);
const ELEMENT_FACE * faces = MeshTopology::GetFaces (HEX); const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (HEX);
for (int j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
{ {
@ -1448,7 +1448,7 @@ namespace netgen
{ {
Element el = mesh[i] ; Element el = mesh[i] ;
HPRefElement & hpel = hpelements[mesh[i].hp_elnr]; HPRefElement & hpel = hpelements[mesh[i].hp_elnr];
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (mesh[i].GetType()); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType());
double dist[3] = {0,0,0}; double dist[3] = {0,0,0};
int ord_dir[3] = {0,0,0}; int ord_dir[3] = {0,0,0};
int edge_dir[12] = {0,0,0,0,0,0,0,0,0,0,0,0}; int edge_dir[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
@ -1520,7 +1520,7 @@ namespace netgen
{ {
Element2d el = mesh[i] ; Element2d el = mesh[i] ;
HPRefElement & hpel = hpelements[mesh[i].hp_elnr]; HPRefElement & hpel = hpelements[mesh[i].hp_elnr];
const ELEMENT_EDGE * edges = MeshTopology::GetEdges (mesh[i].GetType()); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType());
double dist[3] = {0,0,0}; double dist[3] = {0,0,0};
int ord_dir[3] = {0,0,0}; int ord_dir[3] = {0,0,0};
int edge_dir[4] = {0,0,0,0} ; int edge_dir[4] = {0,0,0,0} ;
@ -1637,7 +1637,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
for (int i = 1; i <= mesh.GetNE(); i++) for (int i = 1; i <= mesh.GetNE(); i++)
{ {
const Element & el = mesh.VolumeElement(i); const Element & el = mesh.VolumeElement(i);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (el.GetType()); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (el.GetType());
int nedges = MeshTopology::GetNEdges (el.GetType()); int nedges = MeshTopology::GetNEdges (el.GetType());
for (int j = 0; j < nedges; j++) for (int j = 0; j < nedges; j++)
for (int k = 0; k < nedges; k++) for (int k = 0; k < nedges; k++)

473
libsrc/meshing/topology.cpp
View File

@ -2,6 +2,7 @@
#include "meshing.hpp" #include "meshing.hpp"
#define DEBUG
namespace netgen namespace netgen
{ {
@ -132,6 +133,8 @@ void MeshTopology :: Update()
vert2segment->AddSave (seg[1], i); vert2segment->AddSave (seg[1], i);
} }
if (buildedges) if (buildedges)
{ {
static int timer1 = NgProfiler::CreateTimer ("topology::buildedges"); static int timer1 = NgProfiler::CreateTimer ("topology::buildedges");
@ -162,22 +165,33 @@ void MeshTopology :: Update()
cnt = 0; cnt = 0;
for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++) for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++)
{ {
/* JS, Oct 2009
int pa[2]; int pa[2];
pa[0] = mesh.mlbetweennodes[i].I1(); pa[0] = mesh.mlbetweennodes[i].I1();
pa[1] = mesh.mlbetweennodes[i].I2(); pa[1] = mesh.mlbetweennodes[i].I2();
if (pa[0] > pa[1]) Swap (pa[0], pa[1]); if (pa[0] > pa[1]) Swap (pa[0], pa[1]);
if (pa[0] > 0) if (pa[0] > 0)
cnt.Elem(pa[0])++; cnt.Elem(pa[0])++;
*/
INDEX_2 parents = mesh.mlbetweennodes[i];
parents.Sort();
if (parents[0] >= PointIndex::BASE) cnt[parents[0]]++;
} }
TABLE<int,PointIndex::BASE> vert2vertcoarse (cnt); TABLE<int,PointIndex::BASE> vert2vertcoarse (cnt);
for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++) for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++)
{ {
/*
int pa[2]; int pa[2];
pa[0] = mesh.mlbetweennodes[i].I1(); pa[0] = mesh.mlbetweennodes[i].I1();
pa[1] = mesh.mlbetweennodes[i].I2(); pa[1] = mesh.mlbetweennodes[i].I2();
if (pa[0] > pa[1]) swap (pa[0], pa[1]); if (pa[0] > pa[1]) swap (pa[0], pa[1]);
if (pa[0] > 0) if (pa[0] > 0)
vert2vertcoarse.AddSave1 (pa[0], pa[1]); vert2vertcoarse.AddSave1 (pa[0], pa[1]);
*/
INDEX_2 parents = mesh.mlbetweennodes[i];
parents.Sort();
if (parents[0] > PointIndex::BASE) vert2vertcoarse.AddSave (parents[0], parents[1]);
} }
@ -188,18 +202,18 @@ void MeshTopology :: Update()
ned = edge2vert.Size(); ned = edge2vert.Size();
Array<INDEX_3> missing; Array<INDEX_3> missing;
for (int i = 1; i <= nv; i++) for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++)
{ {
for (int j = 1; j <= vert2edge.EntrySize(i); j++) for (int j = 0; j < vert2edge[i].Size(); j++)
{ {
int ednr = vert2edge.Get(i,j); int ednr = vert2edge[i][j];
int i2 = edge2vert.Get(ednr)[1]; int i2 = edge2vert.Get(ednr)[1];
edgeflag[i2] = i; edgeflag[i2] = i;
edgenr[i2] = ednr; edgenr[i2] = ednr;
} }
for (int j = 1; j <= vert2vertcoarse.EntrySize(i); j++) for (int j = 0; j << vert2vertcoarse[i].Size(); j++)
{ {
int v2 = vert2vertcoarse.Get(i,j); int v2 = vert2vertcoarse[i][j];
if (edgeflag[v2] < i) if (edgeflag[v2] < i)
{ {
ned++; ned++;
@ -209,18 +223,17 @@ void MeshTopology :: Update()
} }
} }
for (int j = 1; j <= vert2element->EntrySize(i); j++) for (int j = 0; j < (*vert2element)[i].Size(); j++)
{ {
int elnr = vert2element->Get(i,j); int elnr = (*vert2element)[i][j];
const Element & el = mesh.VolumeElement (elnr); const Element & el = mesh.VolumeElement (elnr);
int neledges = GetNEdges (el.GetType()); int neledges = GetNEdges (el.GetType());
const ELEMENT_EDGE * eledges = GetEdges (el.GetType()); const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
for (int k = 0; k < neledges; k++) for (int k = 0; k < neledges; k++)
{ {
INDEX_2 edge(el.PNum(eledges[k][0]), INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
el.PNum(eledges[k][1]));
int edgedir = (edge.I1() > edge.I2()); int edgedir = (edge.I1() > edge.I2());
if (edgedir) swap (edge.I1(), edge.I2()); if (edgedir) swap (edge.I1(), edge.I2());
@ -241,18 +254,17 @@ void MeshTopology :: Update()
} }
} }
for (int j = 1; j <= vert2surfelement->EntrySize(i); j++) for (int j = 0; j < (*vert2surfelement)[i].Size(); j++)
{ {
int elnr = vert2surfelement->Get(i,j); int elnr = (*vert2surfelement)[i][j];
const Element2d & el = mesh.SurfaceElement (elnr); const Element2d & el = mesh.SurfaceElement (elnr);
int neledges = GetNEdges (el.GetType()); int neledges = GetNEdges (el.GetType());
const ELEMENT_EDGE * eledges = GetEdges (el.GetType()); const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
for (int k = 0; k < neledges; k++) for (int k = 0; k < neledges; k++)
{ {
INDEX_2 edge(el.PNum(eledges[k][0]), INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
el.PNum(eledges[k][1]));
int edgedir = (edge.I1() > edge.I2()); int edgedir = (edge.I1() > edge.I2());
if (edgedir) swap (edge.I1(), edge.I2()); if (edgedir) swap (edge.I1(), edge.I2());
@ -273,9 +285,9 @@ void MeshTopology :: Update()
} }
} }
for (int j = 1; j <= vert2segment->EntrySize(i); j++) for (int j = 0; j < (*vert2segment)[i].Size(); j++)
{ {
int elnr = vert2segment->Get(i,j); int elnr = (*vert2segment)[i][j];
const Segment & el = mesh.LineSegment (elnr); const Segment & el = mesh.LineSegment (elnr);
INDEX_2 edge(el[0], el[1]); INDEX_2 edge(el[0], el[1]);
@ -306,12 +318,11 @@ void MeshTopology :: Update()
const Element & el = mesh.VolumeElement (i); const Element & el = mesh.VolumeElement (i);
int neledges = GetNEdges (el.GetType()); int neledges = GetNEdges (el.GetType());
const ELEMENT_EDGE * eledges = GetEdges (el.GetType()); const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
for (int k = 0; k < neledges; k++) for (int k = 0; k < neledges; k++)
{ {
INDEX_2 edge(el.PNum(eledges[k][0]), INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
el.PNum(eledges[k][1]));
int edgedir = (edge.I1() > edge.I2()); int edgedir = (edge.I1() > edge.I2());
if (edgedir) swap (edge.I1(), edge.I2()); if (edgedir) swap (edge.I1(), edge.I2());
@ -327,12 +338,11 @@ void MeshTopology :: Update()
const Element2d & el = mesh.SurfaceElement (i); const Element2d & el = mesh.SurfaceElement (i);
int neledges = GetNEdges (el.GetType()); int neledges = GetNEdges (el.GetType());
const ELEMENT_EDGE * eledges = GetEdges (el.GetType()); const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
for (int k = 0; k < neledges; k++) for (int k = 0; k < neledges; k++)
{ {
INDEX_2 edge(el.PNum(eledges[k][0]), INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
el.PNum(eledges[k][1]));
int edgedir = (edge.I1() > edge.I2()); int edgedir = (edge.I1() > edge.I2());
if (edgedir) swap (edge.I1(), edge.I2()); if (edgedir) swap (edge.I1(), edge.I2());
@ -388,13 +398,18 @@ void MeshTopology :: Update()
// generate faces // generate faces
if (buildfaces) // && mesh.GetDimension() == 3) if (buildfaces) // && mesh.GetDimension() == 3)
{ {
int i, j;
static int timer2 = NgProfiler::CreateTimer ("topology::buildfaces"); static int timer2 = NgProfiler::CreateTimer ("topology::buildfaces");
NgProfiler::RegionTimer reg2 (timer2); NgProfiler::RegionTimer reg2 (timer2);
PrintMessage (5, "Update faces "); PrintMessage (5, "Update faces ");
/*
faces.SetSize(ne); faces.SetSize(ne);
surffaces.SetSize(nse); surffaces.SetSize(nse);
@ -403,7 +418,7 @@ void MeshTopology :: Update()
// INDEX_3_HASHTABLE<int> vert2face(ne+nse+1); // INDEX_3_HASHTABLE<int> vert2face(ne+nse+1);
INDEX_3_CLOSED_HASHTABLE<int> vert2face(8*ne+2*nse+nfa+2); INDEX_3_CLOSED_HASHTABLE<int> vert2face(8*ne+2*nse+nfa+2);
for (i = 1; i <= face2vert.Size(); i++) for (int i = 1; i <= face2vert.Size(); i++)
{ {
INDEX_3 f; INDEX_3 f;
f.I1() = face2vert.Get(i)[0]; f.I1() = face2vert.Get(i)[0];
@ -412,16 +427,16 @@ void MeshTopology :: Update()
vert2face.Set (f, i); vert2face.Set (f, i);
} }
for (i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
{ {
const Element & el = mesh.VolumeElement (i); const Element & el = mesh.VolumeElement (i);
int nelfaces = GetNFaces (el.GetType()); int nelfaces = GetNFaces (el.GetType());
const ELEMENT_FACE * elfaces = GetFaces (el.GetType()); const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType());
for (j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
faces.Elem(i)[j] = 0; faces.Elem(i)[j] = 0;
for (j = 0; j < nelfaces; j++) for (int j = 0; j < nelfaces; j++)
if (elfaces[j][3] == 0) if (elfaces[j][3] == 0)
{ // triangle { // triangle
@ -522,14 +537,14 @@ void MeshTopology :: Update()
} }
face2surfel.SetSize(nfa+nse); face2surfel.SetSize(nfa+nse);
for (i = 1; i <= face2surfel.Size(); i++) for (int i = 1; i <= face2surfel.Size(); i++)
face2surfel.Elem(i) = 0; face2surfel.Elem(i) = 0;
for (i = 1; i <= nse; i++) for (int i = 1; i <= nse; i++)
{ {
const Element2d & el = mesh.SurfaceElement (i); const Element2d & el = mesh.SurfaceElement (i);
const ELEMENT_FACE * elfaces = GetFaces (el.GetType()); const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType());
if (elfaces[0][3] == 0) if (elfaces[0][3] == 0)
@ -622,12 +637,6 @@ void MeshTopology :: Update()
// face2vert.SetSize(face2vert.Size()+1); // face2vert.SetSize(face2vert.Size()+1);
INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I3()); INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I3());
face2vert.Append (hface); face2vert.Append (hface);
/*
face2vert.Last()[0] = face4.I1();
face2vert.Last()[1] = face4.I2();
face2vert.Last()[2] = face4.I3();
face2vert.Last()[3] = face4.I3();
*/
} }
surffaces.Elem(i) = 8*(facenum-1)+facedir+1; surffaces.Elem(i) = 8*(facenum-1)+facedir+1;
@ -637,13 +646,13 @@ void MeshTopology :: Update()
surf2volelement.SetSize (nse); surf2volelement.SetSize (nse);
for (i = 1; i <= nse; i++) for (int i = 1; i <= nse; i++)
{ {
surf2volelement.Elem(i)[0] = 0; surf2volelement.Elem(i)[0] = 0;
surf2volelement.Elem(i)[1] = 0; surf2volelement.Elem(i)[1] = 0;
} }
for (i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
for (j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
{ {
int fnum = (faces.Get(i)[j]+7) / 8; int fnum = (faces.Get(i)[j]+7) / 8;
if (fnum > 0 && face2surfel.Elem(fnum)) if (fnum > 0 && face2surfel.Elem(fnum))
@ -656,6 +665,339 @@ void MeshTopology :: Update()
} }
face2vert.SetAllocSize (face2vert.Size()); face2vert.SetAllocSize (face2vert.Size());
*/
faces.SetSize(ne);
surffaces.SetSize(nse);
int oldnfa = face2vert.Size();
cnt = 0;
for (int i = 0; i < face2vert.Size(); i++)
cnt[face2vert[i][0]]++;
TABLE<int,PointIndex::BASE> vert2oldface(cnt);
for (int i = 0; i < face2vert.Size(); i++)
vert2oldface.AddSave (face2vert[i][0], i);
for (int elnr = 0; elnr < ne; elnr++)
for (int j = 0; j < 6; j++)
faces[elnr][j] = 0;
int max_face_on_vertex = 0;
for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++)
{
int onv = vert2oldface[i].Size() + (*vert2element)[i].Size() + (*vert2surfelement)[i].Size();
max_face_on_vertex = max (onv, max_face_on_vertex);
}
for (int pass = 1; pass <= 2; pass++)
{
nfa = oldnfa;
for (int v = PointIndex::BASE; v < nv+PointIndex::BASE; v++)
{
INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10);
for (int j = 0; j < vert2oldface[v].Size(); j++)
{
int fnr = vert2oldface[v][j];
INDEX_3 face (face2vert[fnr].I1(),
face2vert[fnr].I2(),
face2vert[fnr].I3());
vert2face.Set (face, fnr+1);
}
// cout << "inherited faces: " << endl << vert2face << endl;
for (int j = 0; j < (*vert2element)[v].Size(); j++)
{
int elnr = (*vert2element)[v][j];
const Element & el = mesh.VolumeElement (elnr);
int nelfaces = GetNFaces (el.GetType());
const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType());
for (int j = 0; j < nelfaces; j++)
if (elfaces[j][3] == 0)
{ // triangle
int facenum;
int facedir;
INDEX_3 face(el.PNum(elfaces[j][0]),
el.PNum(elfaces[j][1]),
el.PNum(elfaces[j][2]));
facedir = 0;
if (face.I1() > face.I2())
{
swap (face.I1(), face.I2());
facedir += 1;
}
if (face.I2() > face.I3())
{
swap (face.I2(), face.I3());
facedir += 2;
}
if (face.I1() > face.I2())
{
swap (face.I1(), face.I2());
facedir += 4;
}
if (face.I1() != v) continue;
if (vert2face.Used (face))
facenum = vert2face.Get(face);
else
{
nfa++;
vert2face.Set (face, nfa);
facenum = nfa;
INDEX_4 hface(face.I1(),face.I2(),face.I3(),0);
if (pass == 2) face2vert.Append (hface);
}
// cout << "face = " << face << " elnr = " << elnr << endl;
faces.Elem(elnr)[j] = 8*(facenum-1)+facedir+1;
}
else
{
// quad
int facenum;
int facedir;
INDEX_4Q face4(el.PNum(elfaces[j][0]),
el.PNum(elfaces[j][1]),
el.PNum(elfaces[j][2]),
el.PNum(elfaces[j][3]));
facedir = 0;
if (min2 (face4.I1(), face4.I2()) >
min2 (face4.I4(), face4.I3()))
{ // z - flip
facedir += 1;
swap (face4.I1(), face4.I4());
swap (face4.I2(), face4.I3());
}
if (min2 (face4.I1(), face4.I4()) >
min2 (face4.I2(), face4.I3()))
{ // x - flip
facedir += 2;
swap (face4.I1(), face4.I2());
swap (face4.I3(), face4.I4());
}
if (face4.I2() > face4.I4())
{ // diagonal flip
facedir += 4;
swap (face4.I2(), face4.I4());
}
INDEX_3 face(face4.I1(), face4.I2(), face4.I3());
if (face.I1() != v) continue;
if (vert2face.Used (face))
{
facenum = vert2face.Get(face);
}
else
{
nfa++;
vert2face.Set (face, nfa);
facenum = nfa;
INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I4());
if (pass == 2) face2vert.Append (hface);
}
faces.Elem(elnr)[j] = 8*(facenum-1)+facedir+1;
}
}
for (int j = 0; j < (*vert2surfelement)[v].Size(); j++)
{
int elnr = (*vert2surfelement)[v][j];
// cout << "surfelnr = " << elnr << endl;
const Element2d & el = mesh.SurfaceElement (elnr);
const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType());
if (elfaces[0][3] == 0)
{ // triangle
int facenum;
int facedir;
INDEX_3 face(el.PNum(elfaces[0][0]),
el.PNum(elfaces[0][1]),
el.PNum(elfaces[0][2]));
// cout << "face = " << face << endl;
facedir = 0;
if (face.I1() > face.I2())
{
swap (face.I1(), face.I2());
facedir += 1;
}
if (face.I2() > face.I3())
{
swap (face.I2(), face.I3());
facedir += 2;
}
if (face.I1() > face.I2())
{
swap (face.I1(), face.I2());
facedir += 4;
}
if (face.I1() != v) continue;
if (vert2face.Used (face))
facenum = vert2face.Get(face);
else
{
nfa++;
vert2face.Set (face, nfa);
facenum = nfa;
INDEX_4 hface(face.I1(),face.I2(),face.I3(),0);
if (pass == 2) face2vert.Append (hface);
}
// cout << "face = " << face << " selnr = " << elnr << endl;
surffaces.Elem(elnr) = 8*(facenum-1)+facedir+1;
// face2surfel.Elem(facenum) = elnr;
}
else
{
// quad
int facenum;
int facedir;
INDEX_4Q face4(el.PNum(elfaces[0][0]),
el.PNum(elfaces[0][1]),
el.PNum(elfaces[0][2]),
el.PNum(elfaces[0][3]));
facedir = 0;
if (min2 (face4.I1(), face4.I2()) >
min2 (face4.I4(), face4.I3()))
{ // z - orientation
facedir += 1;
swap (face4.I1(), face4.I4());
swap (face4.I2(), face4.I3());
}
if (min2 (face4.I1(), face4.I4()) >
min2 (face4.I2(), face4.I3()))
{ // x - orientation
facedir += 2;
swap (face4.I1(), face4.I2());
swap (face4.I3(), face4.I4());
}
if (face4.I2() > face4.I4())
{
facedir += 4;
swap (face4.I2(), face4.I4());
}
INDEX_3 face(face4.I1(), face4.I2(), face4.I3());
if (face.I1() != v) continue;
if (vert2face.Used (face))
facenum = vert2face.Get(face);
else
{
nfa++;
vert2face.Set (face, nfa);
facenum = nfa;
INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I3());
if (pass == 2) face2vert.Append (hface);
}
surffaces.Elem(elnr) = 8*(facenum-1)+facedir+1;
// face2surfel.Elem(facenum) = elnr;
}
}
}
face2vert.SetAllocSize (nfa);
}
face2surfel.SetSize (nfa);
face2surfel = 0;
for (int i = 1; i <= nse; i++)
face2surfel.Elem(GetSurfaceElementFace(i)) = i;
/*
cout << "build table complete" << endl;
cout << "faces = " << endl;
cout << "face2vert = " << endl << face2vert << endl;
cout << "surffaces = " << endl << surffaces << endl;
cout << "face2surfel = " << endl << face2surfel << endl;
*/
surf2volelement.SetSize (nse);
for (int i = 1; i <= nse; i++)
{
surf2volelement.Elem(i)[0] = 0;
surf2volelement.Elem(i)[1] = 0;
}
for (int i = 1; i <= ne; i++)
for (int j = 0; j < 6; j++)
{
int fnum = (faces.Get(i)[j]+7) / 8;
if (fnum > 0 && face2surfel.Elem(fnum))
{
int sel = face2surfel.Elem(fnum);
surf2volelement.Elem(sel)[1] =
surf2volelement.Elem(sel)[0];
surf2volelement.Elem(sel)[0] = i;
}
}
face2vert.SetAllocSize (face2vert.Size());
// face table complete
/* /*
*testout << "face2vert: "; *testout << "face2vert: ";
@ -672,23 +1014,24 @@ void MeshTopology :: Update()
paralleltop.Reset (); paralleltop.Reset ();
#endif #endif
Array<char> face_els(nfa), face_surfels(nfa); Array<short int> face_els(nfa), face_surfels(nfa);
face_els = 0; face_els = 0;
face_surfels = 0; face_surfels = 0;
Array<int> hfaces; Array<int> hfaces;
for (i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
{ {
GetElementFaces (i, hfaces); GetElementFaces (i, hfaces);
for (j = 0; j < hfaces.Size(); j++) for (int j = 0; j < hfaces.Size(); j++)
face_els[hfaces[j]-1]++; face_els[hfaces[j]-1]++;
} }
for (i = 1; i <= nse; i++) for (int i = 1; i <= nse; i++)
face_surfels[GetSurfaceElementFace (i)-1]++; face_surfels[GetSurfaceElementFace (i)-1]++;
if (ne) if (ne)
{ {
int cnt_err = 0; int cnt_err = 0;
for (i = 0; i < nfa; i++) for (int i = 0; i < nfa; i++)
{ {
/* /*
(*testout) << "face " << i << " has " << int(face_els[i]) << " els, " (*testout) << "face " << i << " has " << int(face_els[i]) << " els, "
@ -746,7 +1089,7 @@ void MeshTopology :: Update()
(*testout) << "illegal face : " << i << endl; (*testout) << "illegal face : " << i << endl;
(*testout) << "points = " << face2vert[i] << endl; (*testout) << "points = " << face2vert[i] << endl;
(*testout) << "pos = "; (*testout) << "pos = ";
for (j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if (face2vert[i].I(j+1) >= 1) if (face2vert[i].I(j+1) >= 1)
(*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " "; (*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " ";
(*testout) << endl; (*testout) << endl;
@ -953,36 +1296,40 @@ void MeshTopology :: GetElementEdges (int elnr, Array<int> & eledges) const
} }
void MeshTopology :: GetElementFaces (int elnr, Array<int> & elfaces, bool withorientation) const void MeshTopology :: GetElementFaces (int elnr, Array<int> & elfaces, bool withorientation) const
{ {
int i;
int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType());
elfaces.SetSize (nfa); elfaces.SetSize (nfa);
for (i = 1; i <= nfa; i++)
if (!withorientation)
for (int i = 1; i <= nfa; i++)
{ {
elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1; elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1;
if(withorientation) }
else
for (int i = 1; i <= nfa; i++)
{ {
elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1;
int orient = (faces.Get(elnr)[i-1]-1) % 8; int orient = (faces.Get(elnr)[i-1]-1) % 8;
if(orient == 1 || orient == 2 || orient == 4 || orient == 7) if(orient == 1 || orient == 2 || orient == 4 || orient == 7)
elfaces.Elem(i) *= -1; elfaces.Elem(i) *= -1;
} }
} }
}
void MeshTopology :: GetElementEdgeOrientations (int elnr, Array<int> & eorient) const void MeshTopology :: GetElementEdgeOrientations (int elnr, Array<int> & eorient) const
{ {
int i;
int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); int ned = GetNEdges (mesh.VolumeElement(elnr).GetType());
eorient.SetSize (ned); eorient.SetSize (ned);
for (i = 1; i <= ned; i++) for (int i = 1; i <= ned; i++)
eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1;
} }
void MeshTopology :: GetElementFaceOrientations (int elnr, Array<int> & forient) const void MeshTopology :: GetElementFaceOrientations (int elnr, Array<int> & forient) const
{ {
int i;
int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType());
forient.SetSize (nfa); forient.SetSize (nfa);
for (i = 1; i <= nfa; i++) for (int i = 1; i <= nfa; i++)
forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8; forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8;
} }
@ -990,14 +1337,13 @@ void MeshTopology :: GetElementFaceOrientations (int elnr, Array<int> & forient)
int MeshTopology :: GetElementEdges (int elnr, int * eledges, int * orient) const int MeshTopology :: GetElementEdges (int elnr, int * eledges, int * orient) const
{ {
int i;
// int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); // int ned = GetNEdges (mesh.VolumeElement(elnr).GetType());
if (mesh.GetDimension()==3 || 1) if (mesh.GetDimension()==3 || 1)
{ {
if (orient) if (orient)
{ {
for (i = 0; i < 12; i++) for (int i = 0; i < 12; i++)
{ {
if (!edges.Get(elnr)[i]) return i; if (!edges.Get(elnr)[i]) return i;
eledges[i] = abs (edges.Get(elnr)[i]); eledges[i] = abs (edges.Get(elnr)[i]);
@ -1006,7 +1352,7 @@ int MeshTopology :: GetElementEdges (int elnr, int * eledges, int * orient) cons
} }
else else
{ {
for (i = 0; i < 12; i++) for (int i = 0; i < 12; i++)
{ {
if (!edges.Get(elnr)[i]) return i; if (!edges.Get(elnr)[i]) return i;
eledges[i] = abs (edges.Get(elnr)[i]); eledges[i] = abs (edges.Get(elnr)[i]);
@ -1041,11 +1387,10 @@ int MeshTopology :: GetElementEdges (int elnr, int * eledges, int * orient) cons
int MeshTopology :: GetElementFaces (int elnr, int * elfaces, int * orient) const int MeshTopology :: GetElementFaces (int elnr, int * elfaces, int * orient) const
{ {
int i;
// int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); // int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType());
if (orient) if (orient)
{ {
for (i = 0; i < 6; i++) for (int i = 0; i < 6; i++)
{ {
if (!faces.Get(elnr)[i]) return i; if (!faces.Get(elnr)[i]) return i;
elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1;
@ -1054,7 +1399,7 @@ int MeshTopology :: GetElementFaces (int elnr, int * elfaces, int * orient) cons
} }
else else
{ {
for (i = 0; i < 6; i++) for (int i = 0; i < 6; i++)
{ {
if (!faces.Get(elnr)[i]) return i; if (!faces.Get(elnr)[i]) return i;
elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1;
@ -1189,7 +1534,7 @@ void MeshTopology :: GetFaceEdges (int fnr, Array<int> & fedges, bool withorient
{ {
const Element & el = mesh.VolumeElement(els[i]); const Element & el = mesh.VolumeElement(els[i]);
int nref_faces = GetNFaces (el.GetType()); int nref_faces = GetNFaces (el.GetType());
const ELEMENT_FACE * ref_faces = GetFaces (el.GetType()); const ELEMENT_FACE * ref_faces = GetFaces1 (el.GetType());
int nfa_ref_edges = GetNEdges (GetFaceType(fnr)); int nfa_ref_edges = GetNEdges (GetFaceType(fnr));
int cntv = 0,fa=-1; int cntv = 0,fa=-1;
@ -1211,7 +1556,7 @@ void MeshTopology :: GetFaceEdges (int fnr, Array<int> & fedges, bool withorient
if(fa>=0) if(fa>=0)
{ {
const ELEMENT_EDGE * fa_ref_edges = GetEdges(GetFaceType(fnr)); const ELEMENT_EDGE * fa_ref_edges = GetEdges1 (GetFaceType(fnr));
fedges.SetSize(nfa_ref_edges); fedges.SetSize(nfa_ref_edges);
GetElementEdges (els[i], eledges); GetElementEdges (els[i], eledges);

199
libsrc/meshing/topology.hpp
View File

@ -62,8 +62,10 @@ public:
static inline int GetNFaces (ELEMENT_TYPE et); static inline int GetNFaces (ELEMENT_TYPE et);
static const Point3d * GetVertices (ELEMENT_TYPE et); static const Point3d * GetVertices (ELEMENT_TYPE et);
inline static const ELEMENT_EDGE * GetEdges (ELEMENT_TYPE et); inline static const ELEMENT_EDGE * GetEdges1 (ELEMENT_TYPE et);
inline static const ELEMENT_FACE * GetFaces (ELEMENT_TYPE et); inline static const ELEMENT_EDGE * GetEdges0 (ELEMENT_TYPE et);
inline static const ELEMENT_FACE * GetFaces1 (ELEMENT_TYPE et);
inline static const ELEMENT_FACE * GetFaces0 (ELEMENT_TYPE et);
int GetSegmentEdge (int segnr) const { return abs(segedges[segnr-1]); } int GetSegmentEdge (int segnr) const { return abs(segedges[segnr-1]); }
@ -295,7 +297,7 @@ int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
const ELEMENT_EDGE * MeshTopology :: GetEdges (ELEMENT_TYPE et) const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
{ {
static int segm_edges[1][2] = static int segm_edges[1][2] =
{ { 1, 2 }}; { { 1, 2 }};
@ -392,7 +394,113 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges (ELEMENT_TYPE et)
} }
const ELEMENT_FACE * MeshTopology :: GetFaces (ELEMENT_TYPE et)
const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
{
static int segm_edges[1][2] =
{ { 0, 1 }};
static int trig_edges[3][2] =
{ { 2, 0 },
{ 1, 2 },
{ 0, 1 }};
static int quad_edges[4][2] =
{ { 0, 1 },
{ 2, 3 },
{ 3, 0 },
{ 1, 2 }};
static int tet_edges[6][2] =
{ { 3, 0 },
{ 3, 1 },
{ 3, 2 },
{ 0, 1 },
{ 0, 2 },
{ 1, 2 }};
static int prism_edges[9][2] =
{ { 2, 0 },
{ 0, 1 },
{ 2, 1 },
{ 5, 3 },
{ 3, 4 },
{ 5, 4 },
{ 2, 5 },
{ 0, 3 },
{ 1, 4 }};
static int pyramid_edges[8][2] =
{ { 0, 1 },
{ 1, 2 },
{ 0, 3 },
{ 3, 2 },
{ 0, 4 },
{ 1, 4 },
{ 2, 4 },
{ 3, 4 }};
static int hex_edges[12][2] =
{
{ 0, 1 },
{ 2, 3 },
{ 3, 0 },
{ 1, 2 },
{ 4, 5 },
{ 6, 7 },
{ 7, 4 },
{ 5, 6 },
{ 0, 4 },
{ 1, 5 },
{ 2, 6 },
{ 3, 7 },
};
switch (et)
{
case SEGMENT:
case SEGMENT3:
return segm_edges;
case TRIG:
case TRIG6:
return trig_edges;
case QUAD:
case QUAD6:
case QUAD8:
return quad_edges;
case TET:
case TET10:
return tet_edges;
case PYRAMID:
return pyramid_edges;
case PRISM:
case PRISM12:
return prism_edges;
case HEX:
return hex_edges;
default:
cerr << "Ng_ME_GetEdges, illegal element type " << et << endl;
}
return 0;
}
const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
{ {
static const int trig_faces[1][4] = static const int trig_faces[1][4] =
{ { 1, 2, 3, 0 } }; { { 1, 2, 3, 0 } };
@ -474,6 +582,89 @@ const ELEMENT_FACE * MeshTopology :: GetFaces (ELEMENT_TYPE et)
const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
{
static const int trig_faces[1][4] =
{ { 0, 1, 2, -1 } };
static const int quad_faces[1][4] =
{ { 0, 1, 2, 3 } };
static const int tet_faces[4][4] =
{ { 3, 1, 2, -1 },
{ 3, 2, 0, -1 },
{ 3, 0, 1, -1 },
{ 0, 2, 1, -1 } };
static const int prism_faces[5][4] =
{
{ 0, 2, 1, -1 },
{ 3, 4, 5, -1 },
{ 2, 0, 3, 5 },
{ 0, 1, 4, 3 },
{ 1, 2, 5, 4 }
};
static const int pyramid_faces[5][4] =
{
{ 0, 1, 4, -1 },
{ 1, 2, 4, -1 },
{ 2, 3, 4, -1 },
{ 3, 0, 4, -1 },
{ 0, 3, 2, 1 }
};
static const int hex_faces[6][4] =
{
{ 0, 3, 2, 1 },
{ 4, 5, 6, 7 },
{ 0, 1, 5, 4 },
{ 1, 2, 6, 5 },
{ 2, 3, 7, 6 },
{ 3, 0, 4, 7 }
};
switch (et)
{
case TRIG:
case TRIG6:
return trig_faces;
case QUAD:
case QUAD6:
case QUAD8:
return quad_faces;
case TET:
case TET10:
return tet_faces;
case PRISM:
case PRISM12:
return prism_faces;
case PYRAMID:
return pyramid_faces;
case SEGMENT:
case SEGMENT3:
case HEX:
return hex_faces;
default:
cerr << "Ng_ME_GetVertices, illegal element type " << et << endl;
}
return 0;
}

8
libsrc/visualization/vsmesh.cpp
View File

@ -1880,7 +1880,7 @@ namespace netgen
if (curv.IsHighOrder()) // && curv.IsElementCurved(ei)) if (curv.IsHighOrder()) // && curv.IsElementCurved(ei))
{ {
const ELEMENT_FACE * faces = MeshTopology :: GetFaces (TET); const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (TET);
const Point3d * vertices = MeshTopology :: GetVertices (TET); const Point3d * vertices = MeshTopology :: GetVertices (TET);
/* /*
@ -2084,7 +2084,7 @@ namespace netgen
CurvedElements & curv = mesh->GetCurvedElements(); CurvedElements & curv = mesh->GetCurvedElements();
if (curv.IsHighOrder()) // && curv.IsElementCurved(ei)) if (curv.IsHighOrder()) // && curv.IsElementCurved(ei))
{ {
const ELEMENT_FACE * faces = MeshTopology :: GetFaces (PRISM); const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (PRISM);
const Point3d * vertices = MeshTopology :: GetVertices (PRISM); const Point3d * vertices = MeshTopology :: GetVertices (PRISM);
Point<3> grid[11][11]; Point<3> grid[11][11];
@ -2469,7 +2469,7 @@ namespace netgen
glEnd (); glEnd ();
*/ */
const ELEMENT_FACE * faces = MeshTopology :: GetFaces (HEX); const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (HEX);
const Point3d * vertices = MeshTopology :: GetVertices (HEX); const Point3d * vertices = MeshTopology :: GetVertices (HEX);
Point<3> grid[11][11]; Point<3> grid[11][11];
@ -2618,7 +2618,7 @@ namespace netgen
if (curv.IsHighOrder()) // && curv.IsElementCurved(ei)) if (curv.IsHighOrder()) // && curv.IsElementCurved(ei))
{ {
const ELEMENT_FACE * faces = MeshTopology :: GetFaces (PYRAMID); const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (PYRAMID);
const Point3d * vertices = MeshTopology :: GetVertices (PYRAMID); const Point3d * vertices = MeshTopology :: GetVertices (PYRAMID);
Point<3> grid[11][11]; Point<3> grid[11][11];