edge/face topology

This commit is contained in:
Joachim Schoeberl 2014-01-07 10:42:39 +00:00
parent 3d8e1d407f
commit ed51959493
7 changed files with 745 additions and 426 deletions

View File

@ -15,11 +15,16 @@
namespace netgen
{
struct T_EDGE2
{
int orient:1;
int nr:31; // 0-based
};
struct T_FACE2
{
int orient:3;
int nr:29; // 0-based
};
class Ng_Element
{
@ -49,20 +54,22 @@ namespace netgen
{
public:
int num;
const int * ptr;
const T_EDGE2 * ptr;
int Size() const { return num; }
int operator[] (int i) const { return abs (ptr[i])-1; }
// int operator[] (int i) const { return abs (ptr[i])-1; }
int operator[] (int i) const { return ptr[i].nr; }
};
class Ng_Faces
{
public:
int num;
const int * ptr;
const T_FACE2 * ptr;
int Size() const { return num; }
int operator[] (int i) const { return (ptr[i]-1) / 8; }
// int operator[] (int i) const { return (ptr[i]-1) / 8; }
int operator[] (int i) const { return ptr[i].nr; }
};
public:
@ -78,10 +85,12 @@ namespace netgen
class Ng_Point
{
public:
double * pt;
public:
Ng_Point (double * apt) : pt(apt) { ; }
double operator[] (int i)
{ return pt[i]; }
operator const double * () { return pt; }
};
@ -152,8 +161,13 @@ namespace netgen
class Mesh * mesh;
public:
Ngx_Mesh(class Mesh * amesh);
Ngx_Mesh () { ; }
Ngx_Mesh(class Mesh * amesh) : mesh(amesh) { ; }
void LoadMesh (const string & filename);
virtual ~Ngx_Mesh();
bool Valid () { return mesh != NULL; }
int GetDimension() const;
int GetNLevels() const;
@ -197,7 +211,7 @@ namespace netgen
template <int DIM>
Ng_Node<DIM> GetNode (int nr);
Ng_Node<DIM> GetNode (int nr) const;
template <int DIM>
@ -208,16 +222,27 @@ namespace netgen
int FindElementOfPoint
(double * p, double * lami,
bool build_searchtrees = false,
int * const indices = NULL, int numind = 0);
int * const indices = NULL, int numind = 0) const;
};
DLL_HEADER Ngx_Mesh * LoadMesh (const string & filename);
}
#ifdef HAVE_NETGEN_SOURCES
#include <meshing.hpp>
namespace netgen
{
#define NGX_INLINE inline
#include <nginterface_v2_impl.hpp>
}
#endif
#endif

View File

@ -0,0 +1,91 @@
NGX_INLINE DLL_HEADER Ng_Point Ngx_Mesh :: GetPoint (int nr) const
{
return Ng_Point (&mesh->Point(nr + PointIndex::BASE)(0));
}
template <>
NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const
{
return (*mesh)[SegmentIndex(nr)].si;
}
template <>
NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<2> (int nr) const
{
int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex();
return mesh->GetFaceDescriptor(ind).BCProperty();
}
template <>
NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<3> (int nr) const
{
return (*mesh)[ElementIndex(nr)].GetIndex();
}
template <>
NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
{
const Segment & el = mesh->LineSegment (SegmentIndex(nr));
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&(el[0]);
ret.vertices.num = 2;
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = 1;
ret.edges.ptr = (T_EDGE2*)mesh->GetTopology().GetSegmentElementEdgesPtr (nr);
ret.faces.num = 0;
ret.faces.ptr = NULL;
return ret;
}
template <>
NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (int nr) const
{
const Element2d & el = mesh->SurfaceElement (SurfaceElementIndex (nr));
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&el[0];
ret.vertices.num = el.GetNV();
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = MeshTopology::GetNEdges (el.GetType());
ret.edges.ptr = (T_EDGE2*)mesh->GetTopology().GetSurfaceElementEdgesPtr (nr);
ret.faces.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetSurfaceElementFacesPtr (nr);
return ret;
}
template <>
NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (int nr) const
{
const Element & el = mesh->VolumeElement (ElementIndex (nr));
Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.points.num = el.GetNP();
ret.points.ptr = (int*)&el[0];
ret.vertices.num = el.GetNV();
ret.vertices.ptr = (int*)&(el[0]);
ret.edges.num = MeshTopology::GetNEdges (el.GetType());
ret.edges.ptr = (T_EDGE2*)mesh->GetTopology().GetElementEdgesPtr (nr);
ret.faces.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetElementFacesPtr (nr);
return ret;
}

View File

@ -18,6 +18,8 @@ namespace netgen
namespace netgen
{
#define NGX_INLINE
#include "nginterface_v2_impl.hpp"
Ngx_Mesh * LoadMesh (const string & filename)
{
@ -26,9 +28,18 @@ namespace netgen
return new Ngx_Mesh (netgen::mesh.Ptr());
}
void Ngx_Mesh :: LoadMesh (const string & filename)
{
netgen::mesh.Ptr() = NULL;
Ng_LoadMesh (filename.c_str());
mesh = netgen::mesh.Ptr();
}
/*
Ngx_Mesh :: Ngx_Mesh (Mesh * amesh)
: mesh(amesh)
{ ; }
*/
Ngx_Mesh :: ~Ngx_Mesh ()
{
@ -71,13 +82,12 @@ namespace netgen
return -1;
}
/*
Ng_Point Ngx_Mesh :: GetPoint (int nr) const
{
Ng_Point ret;
ret.pt = &mesh->Point(nr + PointIndex::BASE)(0);
return ret;
return Ng_Point (&mesh->Point(nr + PointIndex::BASE)(0));
}
*/
template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (int nr) const
{
@ -86,6 +96,7 @@ namespace netgen
return ret;
}
/*
template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
{
const Segment & el = mesh->LineSegment (SegmentIndex(nr));
@ -149,7 +160,7 @@ namespace netgen
return ret;
}
*/
template <>
DLL_HEADER int Ngx_Mesh :: GetElementIndex<0> (int nr) const
@ -157,6 +168,7 @@ namespace netgen
return 0;
}
/*
template <>
DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const
{
@ -175,7 +187,7 @@ namespace netgen
{
return (*mesh)[ElementIndex(nr)].GetIndex();
}
*/
@ -521,14 +533,14 @@ namespace netgen
return mesh->GetTopology().GetNFaces();
}
template <> DLL_HEADER Ng_Node<1> Ngx_Mesh :: GetNode<1> (int nr)
template <> DLL_HEADER Ng_Node<1> Ngx_Mesh :: GetNode<1> (int nr) const
{
Ng_Node<1> node;
node.vertices.ptr = mesh->GetTopology().GetEdgeVerticesPtr(nr);
return node;
}
template <> DLL_HEADER Ng_Node<2> Ngx_Mesh :: GetNode<2> (int nr)
template <> DLL_HEADER Ng_Node<2> Ngx_Mesh :: GetNode<2> (int nr) const
{
Ng_Node<2> node;
node.vertices.ptr = mesh->GetTopology().GetFaceVerticesPtr(nr);
@ -538,10 +550,10 @@ namespace netgen
template <>
DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <2>
DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <2>
(double * p, double * lami,
bool build_searchtree,
int * const indices, int numind)
int * const indices, int numind) const
{
Array<int> dummy(numind);
@ -573,7 +585,7 @@ namespace netgen
DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <3>
(double * p, double * lami,
bool build_searchtree,
int * const indices, int numind)
int * const indices, int numind) const
{
Array<int> dummy(numind);

View File

@ -239,6 +239,24 @@ namespace netgen
values[i+1] = p1;
}
}
template <class S, class T>
void EvaluateScaled (int n, S x, S y, T * values)
{
S p1 = 1.0, p2 = 0.0, p3;
if (n >= 0)
p2 = values[0] = 1.0;
if (n >= 1)
p1 = values[1] = a[0]*y+b[0]*x;
for (int i = 1; i < n; i++)
{
p3 = p2; p2=p1;
p1 = (a[i]*y+b[i]*x)*p2-c[i]*y*y*p3;
values[i+1] = p1;
}
}
};
class JacobiRecPol : public RecPol
@ -331,8 +349,18 @@ namespace netgen
if (n < 3) return;
Tx hx[50], hy[50*50];
ScaledJacobiPolynomial (n-3, x, 1-y, 2, 2, hx);
// ScaledJacobiPolynomial (n-3, x, 1-y, 2, 2, hx);
/*
cout << "scaled jacobi, old: " << endl;
for (int i = 0; i <= n-3; i++)
cout << i << ": " << hx[i] << endl;
*/
jacpols2[2] -> EvaluateScaled (n-3, x, 1-y, hx);
/*
cout << "scaled jacobi, new: " << endl;
for (int i = 0; i <= n-3; i++)
cout << i << ": " << hx[i] << endl;
*/
// for (int ix = 0; ix <= n-3; ix++)
// JacobiPolynomial (n-3, 2*y-1, 2*ix+5, 2, hy+50*ix);
@ -351,7 +379,6 @@ namespace netgen
}
static void CalcTrigShapeDxDy (int n, double x, double y, double * dshape)
{
if (n < 3) return;
@ -366,35 +393,9 @@ namespace netgen
dshape[2*i] = res[i].DValue(0);
dshape[2*i+1] = res[i].DValue(1);
}
/*
if (n < 3) return;
int ndof = (n-1)*(n-2)/2;
double h1[1000], h2[1000];
double eps = 1e-4;
CalcTrigShape (n, x+eps, y, h1);
CalcTrigShape (n, x-eps, y, h2);
for (int i = 0; i < ndof; i++)
dshape[2*i] = (h1[i]-h2[i])/(2*eps);
CalcTrigShape (n, x, y+eps, h1);
CalcTrigShape (n, x, y-eps, h2);
for (int i = 0; i < ndof; i++)
dshape[2*i+1] = (h1[i]-h2[i])/(2*eps);
*/
}
// compute face bubbles up to order n, 0 < y, y-x < 1, x+y < 1
template <class Tx, class Ty, class Tt, class Tr>
static void CalcScaledTrigShape (int n, Tx x, Ty y, Tt t, Tr * shape)
@ -407,7 +408,8 @@ namespace netgen
// ScaledLegendrePolynomial (n-3, (2*y-1), t, hy);
for (int ix = 0; ix <= n-3; ix++)
ScaledJacobiPolynomial (n-3, 2*y-1, t, 2*ix+5, 2, hy+50*ix);
jacpols2[2*ix+5] -> EvaluateScaled (n-3, 2*y-1, t, hy+50*ix);
// ScaledJacobiPolynomial (n-3, 2*y-1, t, 2*ix+5, 2, hy+50*ix);
int ii = 0;
@ -434,44 +436,9 @@ namespace netgen
dshape[3*i+1] = res[i].DValue(1);
dshape[3*i+2] = res[i].DValue(2);
}
/*
double dshape1[6000];
if (n < 3) return;
double hvl[1000], hvr[1000];
int nd = (n-1)*(n-2)/2;
double eps = 1e-6;
CalcScaledTrigShape (n, x-eps, y, t, hvl);
CalcScaledTrigShape (n, x+eps, y, t, hvr);
for (int i = 0; i < nd; i++)
dshape[3*i] = (hvr[i]-hvl[i])/(2*eps);
CalcScaledTrigShape (n, x, y-eps, t, hvl);
CalcScaledTrigShape (n, x, y+eps, t, hvr);
for (int i = 0; i < nd; i++)
dshape[3*i+1] = (hvr[i]-hvl[i])/(2*eps);
CalcScaledTrigShape (n, x, y, t-eps, hvl);
CalcScaledTrigShape (n, x, y, t+eps, hvr);
for (int i = 0; i < nd; i++)
dshape[3*i+2] = (hvr[i]-hvl[i])/(2*eps);
*/
/*
for (int i = 0; i < 3*nd; i++)
if (fabs (dshape[i]-dshape1[i]) > 1e-8 * fabs(dshape[i]) + 1e-6)
{
cerr
cerr << "big difference: " << dshape1[i] << " != " << dshape[i] << endl;
}
*/
}
CurvedElements :: CurvedElements (const Mesh & amesh)
: mesh (amesh)
@ -2168,6 +2135,8 @@ namespace netgen
bool CurvedElements :: IsElementCurved (ElementIndex elnr) const
{
if (mesh[elnr].GetType() != TET) return true;
if (mesh.coarsemesh)
{
const HPRefElement & hpref_el =
@ -2178,7 +2147,7 @@ namespace netgen
const Element & el = mesh[elnr];
ELEMENT_TYPE type = el.GetType();
int nfaces = MeshTopology::GetNFaces (type);
if (nfaces > 4)
{ // not a tet
@ -2225,6 +2194,42 @@ namespace netgen
}
bool CurvedElements :: IsElementHighOrder (ElementIndex elnr) const
{
if (mesh.coarsemesh)
{
const HPRefElement & hpref_el =
(*mesh.hpelements) [mesh[elnr].hp_elnr];
return mesh.coarsemesh->GetCurvedElements().IsElementHighOrder (hpref_el.coarse_elnr);
}
const Element & el = mesh[elnr];
ELEMENT_TYPE type = el.GetType();
ElementInfo info;
info.elnr = elnr;
info.order = order;
info.ndof = info.nv = MeshTopology::GetNPoints (type);
if (info.order > 1)
{
const MeshTopology & top = mesh.GetTopology();
info.nedges = top.GetElementEdges (elnr+1, info.edgenrs, 0);
for (int i = 0; i < info.nedges; i++) info.edgenrs[i]--;
info.nfaces = top.GetElementFaces (elnr+1, info.facenrs, 0);
for (int i = 0; i < info.nfaces; i++) info.facenrs[i]--;
for (int i = 0; i < info.nedges; i++)
if (edgecoeffsindex[info.edgenrs[i]+1] > edgecoeffsindex[info.edgenrs[i]]) return true;
for (int i = 0; i < info.nfaces; i++)
if (facecoeffsindex[info.facenrs[i]+1] > facecoeffsindex[info.facenrs[i]]) return true;
}
return false;
}
@ -3589,6 +3594,8 @@ namespace netgen
double * x, size_t sx,
double * dxdxi, size_t sdxdxi)
{
// static int timer = NgProfiler::CreateTimer ("calcmultipointtrafo, calcshape");
if (mesh.coarsemesh)
{
const HPRefElement & hpref_el =
@ -3675,6 +3682,7 @@ namespace netgen
const Element & el = mesh[elnr];
ELEMENT_TYPE type = el.GetType();
ElementInfo info;
info.elnr = elnr;
info.order = order;

View File

@ -52,6 +52,7 @@ public:
bool IsSegmentCurved (SegmentIndex segnr) const;
bool IsSurfaceElementCurved (SurfaceElementIndex sei) const;
bool IsElementCurved (ElementIndex ei) const;
bool IsElementHighOrder (ElementIndex ei) const;
void CalcSegmentTransformation (double xi, SegmentIndex segnr,

File diff suppressed because it is too large Load Diff

View File

@ -13,6 +13,19 @@
*/
struct T_EDGE
{
int orient:1;
int nr:31; // 0-based
};
struct T_FACE
{
int forient:3;
int fnr:29; // 0-based
};
class MeshTopology
{
const Mesh & mesh;
@ -21,11 +34,11 @@ class MeshTopology
Array<INDEX_2> edge2vert;
Array<INDEX_4> face2vert;
Array<int[12]> edges;
Array<int[6]> faces;
Array<int[4]> surfedges;
Array<int> segedges;
Array<int> surffaces;
Array<T_EDGE[12]> edges;
Array<T_FACE[6]> faces;
Array<T_EDGE[4]> surfedges;
Array<T_EDGE> segedges;
Array<T_FACE> surffaces;
Array<INDEX_2> surf2volelement;
Array<int> face2surfel;
TABLE<ElementIndex,PointIndex::BASE> *vert2element;
@ -66,14 +79,20 @@ public:
inline static const ELEMENT_FACE * GetFaces0 (ELEMENT_TYPE et);
int GetSegmentEdge (int segnr) const { return abs(segedges[segnr-1]); }
int GetSegmentEdgeOrientation (int segnr) const { return sgn(segedges[segnr-1]); }
int GetEdge (SegmentIndex segnr) const { return abs(segedges[segnr])-1; }
// int GetSegmentEdge (int segnr) const { return abs(segedges[segnr-1]); }
// int GetSegmentEdgeOrientation (int segnr) const { return sgn(segedges[segnr-1]); }
// int GetEdge (SegmentIndex segnr) const { return abs(segedges[segnr])-1; }
int GetSegmentEdge (int segnr) const { return segedges[segnr-1].nr+1; }
int GetEdge (SegmentIndex segnr) const { return segedges[segnr].nr; }
void GetSegmentEdge (int segnr, int & enr, int & orient) const
{
/*
enr = abs(segedges.Get(segnr));
orient = segedges.Get(segnr) > 0 ? 1 : -1;
*/
enr = segedges.Get(segnr).nr+1;
orient = segedges.Get(segnr).orient;
}
void GetElementEdges (int elnr, Array<int> & edges) const;
@ -103,12 +122,12 @@ public:
int GetSurfaceElementEdges (int elnr, int * edges, int * orient) const;
const int * GetElementEdgesPtr (int elnr) const { return &edges[elnr][0]; }
const int * GetSurfaceElementEdgesPtr (int selnr) const { return &surfedges[selnr][0]; }
const int * GetSegmentElementEdgesPtr (int selnr) const { return &segedges[selnr]; }
const T_EDGE * GetElementEdgesPtr (int elnr) const { return &edges[elnr][0]; }
const T_EDGE * GetSurfaceElementEdgesPtr (int selnr) const { return &surfedges[selnr][0]; }
const T_EDGE * GetSegmentElementEdgesPtr (int selnr) const { return &segedges[selnr]; }
const int * GetElementFacesPtr (int elnr) const { return &faces[elnr][0]; }
const int * GetSurfaceElementFacesPtr (int selnr) const { return &surffaces[selnr]; }
const T_FACE * GetElementFacesPtr (int elnr) const { return &faces[elnr][0]; }
const T_FACE * GetSurfaceElementFacesPtr (int selnr) const { return &surffaces[selnr]; }
void GetSurface2VolumeElement (int selnr, int & elnr1, int & elnr2) const
@ -141,7 +160,7 @@ public:
int MeshTopology :: GetNVertices (ELEMENT_TYPE et)
inline int MeshTopology :: GetNVertices (ELEMENT_TYPE et)
{
switch (et)
{
@ -172,14 +191,14 @@ int MeshTopology :: GetNVertices (ELEMENT_TYPE et)
case HEX:
return 8;
default:
cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl;
// default:
// cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl;
}
return 0;
}
int MeshTopology :: GetNPoints (ELEMENT_TYPE et)
inline int MeshTopology :: GetNPoints (ELEMENT_TYPE et)
{
switch (et)
{
@ -213,15 +232,15 @@ int MeshTopology :: GetNPoints (ELEMENT_TYPE et)
case HEX:
return 8;
default:
cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl;
// default:
// cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl;
}
return 0;
}
int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
inline int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
{
switch (et)
{
@ -252,14 +271,14 @@ int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
case HEX:
return 12;
default:
cerr << "Ng_ME_GetNEdges, illegal element type " << et << endl;
// default:
// cerr << "Ng_ME_GetNEdges, illegal element type " << et << endl;
}
return 0;
}
int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
inline int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
{
switch (et)
{
@ -290,8 +309,8 @@ int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
case HEX:
return 6;
default:
cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl;
// default:
// cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl;
}
return 0;
}
@ -391,8 +410,8 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
case HEX:
return hex_edges;
default:
cerr << "Ng_ME_GetEdges, illegal element type " << et << endl;
// default:
// cerr << "Ng_ME_GetEdges, illegal element type " << et << endl;
}
return 0;
}
@ -489,8 +508,8 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
case HEX:
return hex_edges;
default:
cerr << "Ng_ME_GetEdges, illegal element type " << et << endl;
// default:
// cerr << "Ng_ME_GetEdges, illegal element type " << et << endl;
}
return 0;
}
@ -504,7 +523,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
{
static const int trig_faces[1][4] =
{ { 1, 2, 3, 0 } };
@ -576,8 +595,8 @@ const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
case HEX:
return hex_faces;
default:
cerr << "Ng_ME_GetVertices, illegal element type " << et << endl;
// default:
// cerr << "Ng_ME_GetVertices, illegal element type " << et << endl;
}
return 0;
}
@ -586,7 +605,7 @@ const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
{
static const int trig_faces[1][4] =
{ { 0, 1, 2, -1 } };
@ -658,8 +677,8 @@ const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
case HEX:
return hex_faces;
default:
cerr << "Ng_ME_GetVertices, illegal element type " << et << endl;
// default:
// cerr << "Ng_ME_GetVertices, illegal element type " << et << endl;
}
return 0;
}