mirror of
https://github.com/NGSolve/netgen.git
synced 2024-09-23 02:02:49 +05:00
365 lines
9.9 KiB
C++
365 lines
9.9 KiB
C++
#include <mystdlib.h>
|
|
|
|
#include "meshing.hpp"
|
|
|
|
namespace netgen
|
|
{
|
|
|
|
AnisotropicClusters :: AnisotropicClusters (const Mesh & amesh)
|
|
: mesh(amesh)
|
|
{
|
|
;
|
|
}
|
|
|
|
AnisotropicClusters :: ~AnisotropicClusters ()
|
|
{
|
|
;
|
|
}
|
|
|
|
void AnisotropicClusters :: Update()
|
|
{
|
|
static Timer timer("clusters");
|
|
// static int timer1 = NgProfiler::CreateTimer ("clusters1");
|
|
// static int timer2 = NgProfiler::CreateTimer ("clusters2");
|
|
// static int timer3 = NgProfiler::CreateTimer ("clusters3");
|
|
RegionTimer reg (timer);
|
|
|
|
const MeshTopology & top = mesh.GetTopology();
|
|
|
|
auto id = this->mesh.GetCommunicator().Rank();
|
|
// auto ntasks = this->mesh.GetCommunicator().Size();
|
|
|
|
bool hasedges = top.HasEdges();
|
|
bool hasfaces = top.HasFaces();
|
|
|
|
if (!hasedges || !hasfaces) return;
|
|
|
|
if (id == 0)
|
|
PrintMessage (3, "Update clusters");
|
|
|
|
nv = mesh.GetNV();
|
|
ned = top.GetNEdges();
|
|
nfa = top.GetNFaces();
|
|
ne = mesh.GetNE();
|
|
int nse = mesh.GetNSE();
|
|
|
|
cluster_reps.SetSize (nv+ned+nfa+ne);
|
|
cluster_reps = -1;
|
|
|
|
NgArray<int> llist (nv+ned+nfa+ne);
|
|
llist = 0;
|
|
|
|
NgArray<int> nnums, ednums, fanums;
|
|
int changed;
|
|
|
|
// NgProfiler::StartTimer(timer1);
|
|
|
|
|
|
/*
|
|
for (int i = 1; i <= ne; i++)
|
|
{
|
|
const Element & el = mesh.VolumeElement(i);
|
|
ELEMENT_TYPE typ = el.GetType();
|
|
|
|
top.GetElementEdges (i, ednums);
|
|
top.GetElementFaces (i, fanums);
|
|
|
|
int elnv = top.GetNVertices (typ);
|
|
int elned = ednums.Size();
|
|
int elnfa = fanums.Size();
|
|
|
|
nnums.SetSize(elnv+elned+elnfa+1);
|
|
for (int j = 1; j <= elnv; j++)
|
|
nnums.Elem(j) = el.PNum(j);
|
|
for (int j = 1; j <= elned; j++)
|
|
nnums.Elem(elnv+j) = nv+ednums.Elem(j);
|
|
for (int j = 1; j <= elnfa; j++)
|
|
nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j);
|
|
nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i;
|
|
|
|
for (int j = 0; j < nnums.Size(); j++)
|
|
cluster_reps.Elem(nnums[j]) = nnums[j];
|
|
}
|
|
*/
|
|
ngcore::ParallelForRange
|
|
(mesh.VolumeElements().Range(),
|
|
[&] (auto myrange)
|
|
{
|
|
NgArray<int> nnums, ednums, fanums;
|
|
for (int i_ : myrange)
|
|
{
|
|
int i = i_+1;
|
|
const Element & el = mesh.VolumeElement(i);
|
|
ELEMENT_TYPE typ = el.GetType();
|
|
|
|
top.GetElementEdges (i, ednums);
|
|
top.GetElementFaces (i, fanums);
|
|
|
|
int elnv = top.GetNVertices (typ);
|
|
int elned = ednums.Size();
|
|
int elnfa = fanums.Size();
|
|
|
|
nnums.SetSize(elnv+elned+elnfa+1);
|
|
for (int j = 1; j <= elnv; j++)
|
|
nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
|
|
for (int j = 1; j <= elned; j++)
|
|
nnums.Elem(elnv+j) = nv+ednums.Elem(j);
|
|
for (int j = 1; j <= elnfa; j++)
|
|
nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j);
|
|
nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i;
|
|
|
|
for (int j = 0; j < nnums.Size(); j++)
|
|
cluster_reps.Elem(nnums[j]) = nnums[j];
|
|
}
|
|
}, ngcore::TasksPerThread(4));
|
|
|
|
// NgProfiler::StopTimer(timer1);
|
|
// NgProfiler::StartTimer(timer2);
|
|
/*
|
|
for (int i = 1; i <= nse; i++)
|
|
{
|
|
const Element2d & el = mesh.SurfaceElement(i);
|
|
ELEMENT_TYPE typ = el.GetType();
|
|
|
|
top.GetSurfaceElementEdges (i, ednums);
|
|
int fanum = top.GetSurfaceElementFace (i);
|
|
|
|
int elnv = top.GetNVertices (typ);
|
|
int elned = ednums.Size();
|
|
|
|
nnums.SetSize(elnv+elned+1);
|
|
for (int j = 1; j <= elnv; j++)
|
|
nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
|
|
for (int j = 1; j <= elned; j++)
|
|
nnums.Elem(elnv+j) = nv+ednums.Elem(j);
|
|
nnums.Elem(elnv+elned+1) = fanum;
|
|
|
|
for (int j = 0; j < nnums.Size(); j++)
|
|
cluster_reps.Elem(nnums[j]) = nnums[j];
|
|
}
|
|
*/
|
|
ngcore::ParallelForRange
|
|
(mesh.SurfaceElements().Range(),
|
|
[&] (auto myrange)
|
|
{
|
|
NgArrayMem<int,9> nnums, ednums;
|
|
for (int i_ : myrange)
|
|
{
|
|
int i = i_+1;
|
|
const Element2d & el = mesh.SurfaceElement(i);
|
|
ELEMENT_TYPE typ = el.GetType();
|
|
|
|
top.GetSurfaceElementEdges (i, ednums);
|
|
int fanum = top.GetSurfaceElementFace (i);
|
|
|
|
int elnv = top.GetNVertices (typ);
|
|
int elned = ednums.Size();
|
|
|
|
nnums.SetSize(elnv+elned+1);
|
|
for (int j = 1; j <= elnv; j++)
|
|
nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
|
|
for (int j = 1; j <= elned; j++)
|
|
nnums.Elem(elnv+j) = nv+ednums.Elem(j);
|
|
nnums.Elem(elnv+elned+1) = fanum;
|
|
|
|
for (int j = 0; j < nnums.Size(); j++)
|
|
cluster_reps.Elem(nnums[j]) = nnums[j];
|
|
}
|
|
}, ngcore::TasksPerThread(4));
|
|
|
|
|
|
// NgProfiler::StopTimer(timer2);
|
|
// NgProfiler::StartTimer(timer3);
|
|
|
|
|
|
static const int hex_cluster[] =
|
|
{
|
|
1, 2, 3, 4, 1, 2, 3, 4,
|
|
5, 6, 7, 8, 5, 6, 7, 8, 1, 2, 3, 4,
|
|
9, 9, 5, 8, 6, 7,
|
|
9
|
|
};
|
|
|
|
static const int prism_cluster[] =
|
|
{
|
|
1, 2, 3, 1, 2, 3,
|
|
4, 5, 6, 4, 5, 6, 3, 1, 2,
|
|
7, 7, 4, 5, 6,
|
|
7
|
|
};
|
|
|
|
static const int pyramid_cluster[] =
|
|
{
|
|
1, 2, 2, 1, 3,
|
|
4, 2, 1, 4, 6, 5, 5, 6,
|
|
7, 5, 7, 6, 4,
|
|
7
|
|
};
|
|
static const int tet_cluster14[] =
|
|
{ 1, 2, 3, 1, 1, 4, 5, 4, 5, 6, 7, 5, 4, 7, 7 };
|
|
|
|
static const int tet_cluster12[] =
|
|
{ 1, 1, 2, 3, 4, 4, 5, 1, 6, 6, 7, 7, 4, 6, 7 };
|
|
|
|
static const int tet_cluster13[] =
|
|
{ 1, 2, 1, 3, 4, 6, 4, 5, 1, 5, 7, 4, 7, 5, 7 };
|
|
|
|
static const int tet_cluster23[] =
|
|
{ 2, 1, 1, 3, 6, 5, 5, 4, 4, 1, 5, 7, 7, 4, 7 };
|
|
|
|
static const int tet_cluster24[] =
|
|
{ 2, 1, 3, 1, 4, 1, 5, 4, 6, 5, 5, 7, 4, 7, 7 };
|
|
|
|
static const int tet_cluster34[] =
|
|
{ 2, 3, 1, 1, 4, 5, 1, 6, 4, 5, 5, 4, 7, 7, 7 };
|
|
|
|
int cnt = 0;
|
|
|
|
do
|
|
{
|
|
static Timer t("update cluster, identify");
|
|
RegionTimer rtr(t);
|
|
cnt++;
|
|
changed = 0;
|
|
|
|
for (int i = 1; i <= ne; i++)
|
|
{
|
|
const Element & el = mesh.VolumeElement(i);
|
|
ELEMENT_TYPE typ = el.GetType();
|
|
|
|
const int * clustertab = NULL;
|
|
switch (typ)
|
|
{
|
|
case PRISM:
|
|
case PRISM12:
|
|
clustertab = prism_cluster;
|
|
break;
|
|
case HEX:
|
|
clustertab = hex_cluster;
|
|
break;
|
|
case PYRAMID:
|
|
clustertab = pyramid_cluster;
|
|
break;
|
|
case TET:
|
|
case TET10:
|
|
if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) ==
|
|
cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE))
|
|
clustertab = tet_cluster12;
|
|
else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) ==
|
|
cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE))
|
|
clustertab = tet_cluster13;
|
|
else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) ==
|
|
cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
|
|
clustertab = tet_cluster14;
|
|
else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) ==
|
|
cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE))
|
|
clustertab = tet_cluster23;
|
|
else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) ==
|
|
cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
|
|
clustertab = tet_cluster24;
|
|
else if (cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE) ==
|
|
cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
|
|
clustertab = tet_cluster34;
|
|
|
|
else
|
|
clustertab = NULL;
|
|
break;
|
|
default:
|
|
clustertab = NULL;
|
|
}
|
|
|
|
if (clustertab)
|
|
{
|
|
top.GetElementEdges (i, ednums);
|
|
top.GetElementFaces (i, fanums);
|
|
|
|
int elnv = top.GetNVertices (typ);
|
|
int elned = ednums.Size();
|
|
int elnfa = fanums.Size();
|
|
|
|
nnums.SetSize(elnv+elned+elnfa+1);
|
|
for (int j = 1; j <= elnv; j++)
|
|
nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
|
|
for (int j = 1; j <= elned; j++)
|
|
nnums.Elem(elnv+j) = nv+ednums.Elem(j);
|
|
for (int j = 1; j <= elnfa; j++)
|
|
nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j);
|
|
nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i;
|
|
|
|
|
|
|
|
for (int j = 0; j < nnums.Size(); j++)
|
|
for (int k = 0; k < j; k++)
|
|
if (clustertab[j] == clustertab[k])
|
|
{
|
|
int jj = nnums[j];
|
|
int kk = nnums[k];
|
|
|
|
if (cluster_reps.Get(kk) < cluster_reps.Get(jj))
|
|
swap (jj,kk);
|
|
|
|
if (cluster_reps.Get(jj) < cluster_reps.Get(kk))
|
|
{
|
|
/*
|
|
cluster_reps.Elem(kk) = cluster_reps.Get(jj);
|
|
changed = 1;
|
|
*/
|
|
|
|
int rep = cluster_reps.Get(jj);
|
|
int next = cluster_reps.Get(kk);
|
|
do
|
|
{
|
|
int cur = next;
|
|
next = llist.Elem(next);
|
|
|
|
cluster_reps.Elem(cur) = rep;
|
|
llist.Elem(cur) = llist.Elem(rep);
|
|
llist.Elem(rep) = cur;
|
|
}
|
|
while (next);
|
|
changed = 1;
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
if (clustertab)
|
|
{
|
|
if (typ == PYRAMID)
|
|
(*testout) << "pyramid";
|
|
else if (typ == PRISM || typ == PRISM12)
|
|
(*testout) << "prism";
|
|
else if (typ == TET || typ == TET10)
|
|
(*testout) << "tet";
|
|
else
|
|
(*testout) << "unknown type" << endl;
|
|
|
|
(*testout) << ", nnums = ";
|
|
for (j = 0; j < nnums.Size(); j++)
|
|
(*testout) << "node " << j << " = " << nnums[j] << ", rep = "
|
|
<< cluster_reps.Get(nnums[j]) << endl;
|
|
}
|
|
*/
|
|
}
|
|
}
|
|
while (changed);
|
|
// NgProfiler::StopTimer(timer3);
|
|
/*
|
|
(*testout) << "cluster reps:" << endl;
|
|
for (i = 1; i <= cluster_reps.Size(); i++)
|
|
{
|
|
(*testout) << i << ": ";
|
|
if (i <= nv)
|
|
(*testout) << "v" << i << " ";
|
|
else if (i <= nv+ned)
|
|
(*testout) << "e" << i-nv << " ";
|
|
else if (i <= nv+ned+nfa)
|
|
(*testout) << "f" << i-nv-ned << " ";
|
|
else
|
|
(*testout) << "c" << i-nv-ned-nfa << " ";
|
|
(*testout) << cluster_reps.Get(i) << endl;
|
|
}
|
|
*/
|
|
}
|
|
}
|