dont't store orientation of edges and faces. This gives a unified memory layout for edges and faces and allows to access facets without copying (in 2D and 3D)

This commit is contained in:
Joachim Schöberl 2016-12-10 19:02:00 +01:00
parent 6513fa3a44
commit e4a6d127fd
4 changed files with 227 additions and 24 deletions

View File

@ -16,13 +16,15 @@ namespace netgen
{
struct T_EDGE2
{
int orient:1;
int nr:31; // 0-based
// int orient:1;
// int nr:31; // 0-based
int nr; // 0-based
};
struct T_FACE2
{
int orient:3;
int nr:29; // 0-based
// int orient:3;
// int nr:29; // 0-based
int nr; // 0-based
};
class Ng_Element
@ -71,6 +73,17 @@ namespace netgen
int operator[] (int i) const { return ptr[i].nr; }
};
class Ng_Facets
{
public:
int num;
const int * ptr;
int Size() const { return num; }
int operator[] (int i) const { return ptr[i]; }
};
public:
NG_ELEMENT_TYPE type;
int index; // material / boundary condition
@ -81,6 +94,7 @@ namespace netgen
Ng_Vertices vertices;
Ng_Edges edges;
Ng_Faces faces;
Ng_Facets facets;
bool is_curved;
};

View File

@ -86,6 +86,20 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
ret.faces.num = 0;
ret.faces.ptr = NULL;
if (mesh->GetDimension() == 2)
{
ret.facets.num = 1;
ret.facets.ptr = (int*)ret.edges.ptr;
}
else
{
ret.facets.num = 0;
ret.facets.ptr = nullptr;
// not working as long as vertices are 1-based
// ret.facets.num = 2;
// ret.facets.ptr = (int*)&(el[0]);
}
// ret.is_curved = mesh->GetCurvedElements().IsSegmentCurved(nr);
ret.is_curved = el.IsCurved();
@ -117,6 +131,16 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (int nr) const
ret.faces.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetSurfaceElementFacesPtr (nr);
if (mesh->GetDimension() == 3)
{
ret.facets.num = ret.faces.num;
ret.facets.ptr = (int*)ret.faces.ptr;
}
else
{
ret.facets.num = ret.edges.num;
ret.facets.ptr = (int*)ret.edges.ptr;
}
ret.is_curved = el.IsCurved();
return ret;
}
@ -142,6 +166,9 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (int nr) const
ret.faces.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetElementFacesPtr (nr);
ret.facets.num = ret.faces.num;
ret.facets.ptr = (int*)ret.faces.ptr;
ret.is_curved = el.IsCurved();
return ret;
}

View File

@ -608,15 +608,15 @@ namespace netgen
{
case 3:
edges[elnr][loc_edge].nr = edgenum;
edges[elnr][loc_edge].orient = edgedir;
// edges[elnr][loc_edge].orient = edgedir;
break;
case 2:
surfedges[elnr][loc_edge].nr = edgenum;
surfedges[elnr][loc_edge].orient = edgedir;
// surfedges[elnr][loc_edge].orient = edgedir;
break;
case 1:
segedges[elnr].nr = edgenum;
segedges[elnr].orient = edgedir;
// segedges[elnr].orient = edgedir;
break;
}
});
@ -783,12 +783,12 @@ namespace netgen
if (volume)
{
faces[elnr][j].fnr = facenum;
faces[elnr][j].forient = facedir;
// faces[elnr][j].forient = facedir;
}
else
{
surffaces[elnr].fnr = facenum;
surffaces[elnr].forient = facedir;
// surffaces[elnr].forient = facedir;
}
});
}
@ -1346,7 +1346,8 @@ namespace netgen
eorient.SetSize (ned);
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].orient) ? -1 : 1;
// eorient.Elem(i) = (edges.Get(elnr)[i-1].orient) ? -1 : 1;
eorient.Elem(i) = GetElementEdgeOrientation (elnr, i-1) ? -1 : 1;
}
void MeshTopology :: GetElementFaceOrientations (int elnr, Array<int> & forient) const
@ -1354,8 +1355,9 @@ namespace netgen
int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType());
forient.SetSize (nfa);
for (int i = 1; i <= nfa; i++)
forient.Elem(i) = faces.Get(elnr)[i-1].forient;
// forient.Elem(i) = faces.Get(elnr)[i-1].forient;
// forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8;
forient.Elem(i) = GetElementFaceOrientation(elnr, i-1);
}
@ -1377,7 +1379,8 @@ namespace netgen
*/
if (edges.Get(elnr)[i].nr == -1) return i;
eledges[i] = edges.Get(elnr)[i].nr+1;
orient[i] = edges.Get(elnr)[i].orient ? -1 : 1;
// orient[i] = edges.Get(elnr)[i].orient ? -1 : 1;
orient[i] = GetElementEdgeOrientation(elnr, i) ? -1 : 1;
}
}
else
@ -1432,7 +1435,8 @@ namespace netgen
*/
if (faces.Get(elnr)[i].fnr == -1) return i;
elfaces[i] = faces.Get(elnr)[i].fnr+1;
orient[i] = faces.Get(elnr)[i].forient;
// orient[i] = faces.Get(elnr)[i].forient;
orient[i] = GetElementFaceOrientation (elnr, i);
}
}
else
@ -1484,13 +1488,15 @@ namespace netgen
eorient.SetSize (ned);
for (int i = 0; i < ned; i++)
// eorient[i] = (surfedges.Get(elnr)[i] > 0) ? 1 : -1;
eorient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1;
// eorient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1;
eorient[i] = GetSurfaceElementEdgeOrientation(elnr, i) ? -1 : 1;
}
int MeshTopology :: GetSurfaceElementFaceOrientation (int elnr) const
{
// return (surffaces.Get(elnr)-1) % 8;
return surffaces.Get(elnr).forient;
// return surffaces.Get(elnr).forient;
return GetSurfaceElementFaceOrientation2(elnr);
}
int MeshTopology :: GetSurfaceElementEdges (int elnr, int * eledges, int * orient) const
@ -1509,7 +1515,8 @@ namespace netgen
*/
if (surfedges.Get(elnr)[i].nr == -1) return i;
eledges[i] = surfedges.Get(elnr)[i].nr+1;
orient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1;
// orient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1;
orient[i] = GetSurfaceElementEdgeOrientation(elnr, i) ? -1 : 1;
}
}
@ -1536,12 +1543,159 @@ namespace netgen
*/
eledges[0] = segedges.Get(elnr).nr+1;
if (orient)
orient[0] = segedges.Get(elnr).orient ? -1 : 1;
// orient[0] = segedges.Get(elnr).orient ? -1 : 1;
orient[0] = GetSegmentEdgeOrientation(elnr) ? -1 : 1;
}
return 1;
}
int MeshTopology :: GetElementEdgeOrientation (int elnr, int locedgenr) const
{
const Element & el = mesh.VolumeElement (elnr);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType());
int k = locedgenr;
INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
int edgedir = (edge.I1() > edge.I2());
return edgedir;
}
int MeshTopology :: GetElementFaceOrientation (int elnr, int locfacenr) const
{
const Element & el = mesh.VolumeElement (elnr);
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType());
int j = locfacenr;
if (elfaces[j][3] < 0)
{ // triangle
INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]],
el[elfaces[j][2]], 0);
int 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; }
return facedir;
}
else
{
// quad
int facenum;
INDEX_4 face4(el[elfaces[j][0]], el[elfaces[j][1]],
el[elfaces[j][2]], el[elfaces[j][3]]);
int 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());
}
return facedir;
}
}
int MeshTopology :: GetSurfaceElementEdgeOrientation (int elnr, int locedgenr) const
{
const Element2d & el = mesh.SurfaceElement (elnr);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType());
int k = locedgenr;
INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
int edgedir = (edge.I1() > edge.I2());
return edgedir;
}
int MeshTopology :: GetSurfaceElementFaceOrientation2 (int elnr) const
{
const Element2d & el = mesh.SurfaceElement (elnr);
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType());
int j = 0;
if (elfaces[j][3] < 0)
{ // triangle
INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]],
el[elfaces[j][2]], 0);
int 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; }
return facedir;
}
else
{
// quad
int facenum;
INDEX_4 face4(el[elfaces[j][0]], el[elfaces[j][1]],
el[elfaces[j][2]], el[elfaces[j][3]]);
int 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());
}
return facedir;
}
}
int MeshTopology :: GetSegmentEdgeOrientation (int elnr) const
{
const Segment & el = mesh.LineSegment (elnr);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType());
int k = 0;
INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
int edgedir = (edge.I1() > edge.I2());
return edgedir;
}
void MeshTopology :: GetFaceVertices (int fnr, Array<int> & vertices) const
{
vertices.SetSize(4);

View File

@ -17,14 +17,14 @@ namespace netgen
struct T_EDGE
{
int orient:1;
int nr:31; // 0-based
// int orient:1;
int nr; // 0-based
};
struct T_FACE
{
int forient:3;
int fnr:29; // 0-based
// int forient:3;
int fnr; // 0-based
};
@ -87,7 +87,8 @@ public:
void GetSegmentEdge (int segnr, int & enr, int & orient) const
{
enr = segedges.Get(segnr).nr+1;
orient = segedges.Get(segnr).orient;
// orient = segedges.Get(segnr).orient;
orient = GetSegmentEdgeOrientation(segnr);
}
void GetElementEdges (int elnr, Array<int> & edges) const;
@ -98,6 +99,13 @@ public:
int GetElementEdges (int elnr, int * edges, int * orient) const;
int GetElementFaces (int elnr, int * faces, int * orient) const;
int GetElementEdgeOrientation (int elnr, int locedgenr) const; // old style
int GetElementFaceOrientation (int elnr, int locfacenr) const; // old style
int GetSurfaceElementEdgeOrientation (int elnr, int locedgenr) const; // old style
int GetSurfaceElementFaceOrientation2 (int elnr) const; // old style
int GetSegmentEdgeOrientation (int elnr) const; // old style
void GetFaceVertices (int fnr, Array<int> & vertices) const;
void GetFaceVertices (int fnr, int * vertices) const;
void GetEdgeVertices (int enr, int & v1, int & v2) const;