mirror of
https://github.com/NGSolve/netgen.git
synced 2024-11-11 16:49:16 +05:00
Merge branch 'master' into stl_refine_fix
This commit is contained in:
Christopher Lackner
commit
9b92e754f2
@ -2,10 +2,271 @@
|
||||
#include "meshing.hpp"
|
||||
|
||||
|
||||
|
||||
namespace netgen
|
||||
{
|
||||
|
||||
template<int dim, typename T=INDEX>
|
||||
class BTree
|
||||
{
|
||||
public:
|
||||
// Number of entries per leaf
|
||||
static constexpr int N = 100;
|
||||
|
||||
struct Node;
|
||||
|
||||
struct Leaf
|
||||
{
|
||||
Point<2*dim> p[N];
|
||||
T index[N];
|
||||
int n_elements;
|
||||
|
||||
Leaf() : n_elements(0) {}
|
||||
|
||||
void Add( Array<Leaf*, INDEX> &leaf_index, const Point<2*dim> &ap, T aindex )
|
||||
{
|
||||
p[n_elements] = ap;
|
||||
index[n_elements] = aindex;
|
||||
n_elements++;
|
||||
if(leaf_index.Size()<=aindex)
|
||||
leaf_index.SetSize(2*aindex);
|
||||
leaf_index[aindex] = this;
|
||||
}
|
||||
};
|
||||
|
||||
struct Node
|
||||
{
|
||||
union
|
||||
{
|
||||
Node *children[2];
|
||||
Leaf *leaf;
|
||||
};
|
||||
double sep;
|
||||
int level;
|
||||
|
||||
Node()
|
||||
: children{nullptr,nullptr}
|
||||
{ }
|
||||
|
||||
~Node()
|
||||
{ }
|
||||
|
||||
Leaf *GetLeaf() const
|
||||
{
|
||||
return children[1] ? nullptr : leaf;
|
||||
}
|
||||
};
|
||||
|
||||
private:
|
||||
Node root;
|
||||
Array<Leaf*, INDEX> leaf_index;
|
||||
Point<dim> global_min, global_max;
|
||||
double tol;
|
||||
size_t n_leaves;
|
||||
size_t n_nodes;
|
||||
BlockAllocator ball_nodes;
|
||||
BlockAllocator ball_leaves;
|
||||
|
||||
public:
|
||||
|
||||
BTree (const Point<dim> & pmin, const Point<dim> & pmax)
|
||||
: global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf))
|
||||
{
|
||||
root.leaf = (Leaf*) ball_leaves.Alloc(); new (root.leaf) Leaf();
|
||||
root.level = 0;
|
||||
tol = 1e-7 * Dist(pmax, pmin);
|
||||
}
|
||||
|
||||
size_t GetNLeaves()
|
||||
{
|
||||
return n_leaves;
|
||||
}
|
||||
|
||||
size_t GetNNodes()
|
||||
{
|
||||
return n_nodes;
|
||||
}
|
||||
|
||||
template<typename TFunc>
|
||||
void GetFirstIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
|
||||
NgArray<T> & pis, TFunc func=[](auto pi){return false;}) const
|
||||
{
|
||||
static Timer timer("BTree::GetIntersecting"); RegionTimer rt(timer);
|
||||
static Timer timer1("BTree::GetIntersecting-LinearSearch");
|
||||
static Array<const Node*> stack(1000);
|
||||
static Array<int> dir_stack(1000);
|
||||
|
||||
pis.SetSize(0);
|
||||
|
||||
Point<2*dim> tpmin, tpmax;
|
||||
|
||||
for (size_t i : IntRange(dim))
|
||||
{
|
||||
tpmin(i) = global_min(i);
|
||||
tpmax(i) = pmax(i)+tol;
|
||||
|
||||
tpmin(i+dim) = pmin(i)-tol;
|
||||
tpmax(i+dim) = global_max(i);
|
||||
}
|
||||
|
||||
stack.SetSize(0);
|
||||
stack.Append(&root);
|
||||
dir_stack.SetSize(0);
|
||||
dir_stack.Append(0);
|
||||
|
||||
while(stack.Size())
|
||||
{
|
||||
const Node *node = stack.Last();
|
||||
stack.DeleteLast();
|
||||
|
||||
int dir = dir_stack.Last();
|
||||
dir_stack.DeleteLast();
|
||||
|
||||
if(Leaf *leaf = node->GetLeaf())
|
||||
{
|
||||
// RegionTimer rt1(timer1);
|
||||
for (auto i : IntRange(leaf->n_elements))
|
||||
{
|
||||
bool intersect = true;
|
||||
const auto p = leaf->p[i];
|
||||
|
||||
for (int d = 0; d < dim; d++)
|
||||
if (p[d] > tpmax[d])
|
||||
intersect = false;
|
||||
for (int d = dim; d < 2*dim; d++)
|
||||
if (p[d] < tpmin[d])
|
||||
intersect = false;
|
||||
if(intersect)
|
||||
{
|
||||
pis.Append (leaf->index[i]);
|
||||
if(func(leaf->index[i])) return;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
int newdir = dir+1;
|
||||
if(newdir==2*dim) newdir = 0;
|
||||
if (tpmin[dir] <= node->sep)
|
||||
{
|
||||
stack.Append(node->children[0]);
|
||||
dir_stack.Append(newdir);
|
||||
}
|
||||
if (tpmax[dir] >= node->sep)
|
||||
{
|
||||
stack.Append(node->children[1]);
|
||||
dir_stack.Append(newdir);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
|
||||
NgArray<T> & pis) const
|
||||
{
|
||||
GetFirstIntersecting(pmin, pmax, pis, [](auto pi){return false;});
|
||||
}
|
||||
|
||||
void Insert (const Point<dim> & pmin, const Point<dim> & pmax, T pi)
|
||||
{
|
||||
// static Timer timer("BTree::Insert"); RegionTimer rt(timer);
|
||||
int dir = 0;
|
||||
Point<2*dim> p;
|
||||
for (auto i : IntRange(dim))
|
||||
{
|
||||
p(i) = pmin[i];
|
||||
p(i+dim) = pmax[i];
|
||||
}
|
||||
|
||||
Node * node = &root;
|
||||
Leaf * leaf = node->GetLeaf();
|
||||
|
||||
// search correct leaf to add point
|
||||
while(!leaf)
|
||||
{
|
||||
node = p[dir] < node->sep ? node->children[0] : node->children[1];
|
||||
dir++;
|
||||
if(dir==2*dim) dir = 0;
|
||||
leaf = node->GetLeaf();
|
||||
}
|
||||
|
||||
// add point to leaf
|
||||
if(leaf->n_elements < N)
|
||||
leaf->Add(leaf_index, p,pi);
|
||||
else // assume leaf->n_elements == N
|
||||
{
|
||||
// add two new nodes and one new leaf
|
||||
int n_elements = leaf->n_elements;
|
||||
ArrayMem<double, N> coords(n_elements);
|
||||
ArrayMem<int, N> order(n_elements);
|
||||
|
||||
// separate points in two halves, first sort all coordinates in direction dir
|
||||
for (auto i : IntRange(n_elements))
|
||||
{
|
||||
order[i] = i;
|
||||
coords[i] = leaf->p[i][dir];
|
||||
}
|
||||
|
||||
QuickSortI(coords, order);
|
||||
int isplit = N/2;
|
||||
Leaf *leaf1 = (Leaf*) ball_leaves.Alloc(); new (leaf1) Leaf();
|
||||
Leaf *leaf2 = (Leaf*) ball_leaves.Alloc(); new (leaf2) Leaf();
|
||||
|
||||
for (auto i : order.Range(isplit))
|
||||
leaf1->Add(leaf_index, leaf->p[i], leaf->index[i] );
|
||||
for (auto i : order.Range(isplit, N))
|
||||
leaf2->Add(leaf_index, leaf->p[i], leaf->index[i] );
|
||||
|
||||
Node *node1 = (Node*) ball_nodes.Alloc(); new (node1) Node();
|
||||
node1->leaf = leaf1;
|
||||
node1->level = node->level+1;
|
||||
|
||||
Node *node2 = (Node*) ball_nodes.Alloc(); new (node2) Node();
|
||||
node2->leaf = leaf2;
|
||||
node2->level = node->level+1;
|
||||
|
||||
node->children[0] = node1;
|
||||
node->children[1] = node2;
|
||||
node->sep = 0.5 * (leaf->p[order[isplit-1]][dir] + leaf->p[order[isplit]][dir]);
|
||||
|
||||
// add new point to one of the new leaves
|
||||
if (p[dir] < node->sep)
|
||||
leaf1->Add( leaf_index, p, pi );
|
||||
else
|
||||
leaf2->Add( leaf_index, p, pi );
|
||||
|
||||
ball_leaves.Free(leaf);
|
||||
n_leaves++;
|
||||
n_nodes+=2;
|
||||
}
|
||||
}
|
||||
|
||||
void DeleteElement (T pi)
|
||||
{
|
||||
// static Timer timer("BTree::DeleteElement"); RegionTimer rt(timer);
|
||||
Leaf *leaf = leaf_index[pi];
|
||||
auto & n_elements = leaf->n_elements;
|
||||
auto & index = leaf->index;
|
||||
auto & p = leaf->p;
|
||||
|
||||
for (auto i : IntRange(n_elements))
|
||||
{
|
||||
if(index[i] == pi)
|
||||
{
|
||||
n_elements--;
|
||||
if(i!=n_elements)
|
||||
{
|
||||
index[i] = index[n_elements];
|
||||
p[i] = p[n_elements];
|
||||
}
|
||||
return;
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
typedef BTree<3> TBoxTree;
|
||||
// typedef BoxTree<3> TBoxTree;
|
||||
|
||||
static const int deltetfaces[][3] =
|
||||
{ { 1, 2, 3 },
|
||||
{ 2, 0, 3 },
|
||||
@ -214,25 +475,41 @@ namespace netgen
|
||||
void AddDelaunayPoint (PointIndex newpi, const Point3d & newp,
|
||||
NgArray<DelaunayTet> & tempels,
|
||||
Mesh & mesh,
|
||||
BoxTree<3> & tettree,
|
||||
TBoxTree & tettree,
|
||||
MeshNB & meshnb,
|
||||
NgArray<Point<3> > & centers, NgArray<double> & radi2,
|
||||
NgArray<int> & connected, NgArray<int> & treesearch,
|
||||
NgArray<int> & freelist, SphereList & list,
|
||||
IndexSet & insphere, IndexSet & closesphere)
|
||||
{
|
||||
// static Timer t("Meshing3::AddDelaunayPoint"); RegionTimer reg(t);
|
||||
static Timer t("Meshing3::AddDelaunayPoint"); RegionTimer reg(t);
|
||||
static Timer tsearch("addpoint, search");
|
||||
static Timer tinsert("addpoint, insert");
|
||||
|
||||
static Timer t0("addpoint, 0");
|
||||
static Timer t1("addpoint, 1");
|
||||
static Timer t2("addpoint, 2");
|
||||
static Timer t3("addpoint, 3");
|
||||
/*
|
||||
find any sphere, such that newp is contained in
|
||||
*/
|
||||
DelaunayTet el;
|
||||
// DelaunayTet el;
|
||||
int cfelind = -1;
|
||||
|
||||
const Point<3> * pp[4];
|
||||
Point<3> pc;
|
||||
Point3d tpmin, tpmax;
|
||||
|
||||
tettree.GetIntersecting (newp, newp, treesearch);
|
||||
|
||||
|
||||
/*
|
||||
// stop search if intersecting point is close enough to center
|
||||
tettree.GetFirstIntersecting (newp, newp, treesearch, [&](const auto pi)
|
||||
{
|
||||
double quot = Dist2 (centers.Get(pi), newp);
|
||||
return (quot < 0.917632 * radi2.Get(pi));
|
||||
} );
|
||||
// tettree.GetIntersecting (newp, newp, treesearch);
|
||||
|
||||
double quot,minquot(1e20);
|
||||
|
||||
@ -249,7 +526,28 @@ namespace netgen
|
||||
break;
|
||||
}
|
||||
}
|
||||
*/
|
||||
|
||||
tsearch.Start();
|
||||
double minquot{1e20};
|
||||
tettree.GetFirstIntersecting
|
||||
(newp, newp, treesearch, [&](const auto pi)
|
||||
{
|
||||
double quot = Dist2 (centers.Get(pi), newp) / radi2.Get(pi);
|
||||
if (quot < 0.917632)
|
||||
{
|
||||
cfelind = pi;
|
||||
return true;
|
||||
}
|
||||
|
||||
if((cfelind == -1 || quot < 0.99*minquot) && quot < 1)
|
||||
{
|
||||
minquot = quot;
|
||||
cfelind = pi;
|
||||
}
|
||||
return false;
|
||||
} );
|
||||
tsearch.Stop();
|
||||
|
||||
if (cfelind == -1)
|
||||
{
|
||||
@ -275,6 +573,7 @@ namespace netgen
|
||||
int changed = 1;
|
||||
int nstarti = 1, starti;
|
||||
|
||||
t0.Start();
|
||||
while (changed)
|
||||
{
|
||||
changed = 0;
|
||||
@ -346,7 +645,7 @@ namespace netgen
|
||||
Vec<3> v1 = p2-p1;
|
||||
Vec<3> v2 = p3-p1;
|
||||
Vec<3> n = Cross (v1, v2);
|
||||
n /= n.Length();
|
||||
n /= n.Length();
|
||||
|
||||
if (n * Vec3d (p1, mesh.Point (tempels.Get(helind)[k])) > 0)
|
||||
n *= -1;
|
||||
@ -363,9 +662,12 @@ namespace netgen
|
||||
}
|
||||
} // while (changed)
|
||||
|
||||
t0.Stop();
|
||||
t1.Start();
|
||||
// NgArray<Element> newels;
|
||||
NgArray<DelaunayTet> newels;
|
||||
|
||||
static NgArray<DelaunayTet> newels;
|
||||
newels.SetSize(0);
|
||||
|
||||
Element2d face(TRIG);
|
||||
|
||||
for (int celind : insphere.GetArray())
|
||||
@ -389,7 +691,7 @@ namespace netgen
|
||||
Vec<3> v2 = mesh[face[2]] - mesh[face[0]];
|
||||
Vec<3> n = Cross (v1, v2);
|
||||
|
||||
n.Normalize();
|
||||
n.Normalize();
|
||||
if (n * Vec3d(mesh.Point (face[0]),
|
||||
mesh.Point (tempels.Get(celind)[k]))
|
||||
> 0)
|
||||
@ -410,6 +712,8 @@ namespace netgen
|
||||
}
|
||||
}
|
||||
|
||||
t1.Stop();
|
||||
t2.Start();
|
||||
meshnb.ResetFaceHT (10*insphere.GetArray().Size()+1);
|
||||
|
||||
for (auto celind : insphere.GetArray())
|
||||
@ -433,7 +737,8 @@ namespace netgen
|
||||
fabs (Dist2 (centers.Get (ind), newp) - radi2.Get(ind)) < 1e-8 )
|
||||
hasclose = true;
|
||||
}
|
||||
|
||||
t2.Stop();
|
||||
t3.Start();
|
||||
for (int j = 1; j <= newels.Size(); j++)
|
||||
{
|
||||
const auto & newel = newels.Get(j);
|
||||
@ -507,8 +812,11 @@ namespace netgen
|
||||
tpmax.SetToMax (*pp[k]);
|
||||
}
|
||||
tpmax = tpmax + 0.01 * (tpmax - tpmin);
|
||||
// tinsert.Start();
|
||||
tettree.Insert (tpmin, tpmax, nelind);
|
||||
// tinsert.Stop();
|
||||
}
|
||||
t3.Stop();
|
||||
}
|
||||
|
||||
|
||||
@ -586,7 +894,7 @@ namespace netgen
|
||||
pmin2 = pmin2 + 0.1 * (pmin2 - pmax2);
|
||||
pmax2 = pmax2 + 0.1 * (pmax2 - pmin2);
|
||||
|
||||
BoxTree<3> tettree(pmin2, pmax2);
|
||||
TBoxTree tettree(pmin2, pmax2);
|
||||
|
||||
|
||||
tempels.Append (startel);
|
||||
@ -619,7 +927,9 @@ namespace netgen
|
||||
|
||||
// "random" reordering of points (speeds a factor 3 - 5 !!!)
|
||||
NgArray<PointIndex, PointIndex::BASE, PointIndex> mixed(np);
|
||||
int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
|
||||
// int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
|
||||
// int prims[] = { 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
|
||||
int prims[] = { 211, 223, 227, 229, 233, 239, 241, 251, 257, 263 };
|
||||
int prim;
|
||||
|
||||
{
|
||||
@ -670,6 +980,8 @@ namespace netgen
|
||||
|
||||
PrintMessage (3, "Points: ", cntp);
|
||||
PrintMessage (3, "Elements: ", tempels.Size());
|
||||
PrintMessage (3, "Tree data entries per element: ", 1.0*tettree.N*tettree.GetNLeaves() / tempels.Size());
|
||||
PrintMessage (3, "Tree nodes per element: ", 1.0*tettree.GetNNodes() / tempels.Size());
|
||||
// (*mycout) << cntp << " / " << tempels.Size() << " points/elements" << endl;
|
||||
|
||||
/*
|
||||
|
||||
@ -891,11 +891,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
|
||||
multithread.task = savetask;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
void MeshOptimize3d :: SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
const NgBitArray * working_elements)
|
||||
{
|
||||
static Timer t("MeshOptimize3d::SwapImprove"); RegionTimer reg(t);
|
||||
@ -1767,6 +1763,806 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
|
||||
multithread.task = savetask;
|
||||
}
|
||||
|
||||
|
||||
bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
const NgBitArray * working_elements,
|
||||
TABLE<ElementIndex,PointIndex::BASE> & elementsonnode,
|
||||
INDEX_3_HASHTABLE<int> & faces,
|
||||
PointIndex pi1, PointIndex pi2, bool check_only)
|
||||
{
|
||||
PointIndex pi3(PointIndex::INVALID), pi4(PointIndex::INVALID),
|
||||
pi5(PointIndex::INVALID), pi6(PointIndex::INVALID);
|
||||
|
||||
double bad1, bad2, bad3;
|
||||
|
||||
Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
|
||||
Element el1(TET), el2(TET), el3(TET), el4(TET);
|
||||
Element el1b(TET), el2b(TET), el3b(TET), el4b(TET);
|
||||
ArrayMem<ElementIndex, 20> hasbothpoints;
|
||||
|
||||
bool do_swap = false;
|
||||
if (pi2 < pi1) Swap (pi1, pi2);
|
||||
|
||||
if (mesh.BoundaryEdge (pi1, pi2)) return false;
|
||||
|
||||
|
||||
hasbothpoints.SetSize (0);
|
||||
for (int k = 0; k < elementsonnode[pi1].Size(); k++)
|
||||
{
|
||||
bool has1 = 0, has2 = 0;
|
||||
ElementIndex elnr = elementsonnode[pi1][k];
|
||||
const Element & elem = mesh[elnr];
|
||||
|
||||
if (elem.IsDeleted()) return false;
|
||||
|
||||
for (int l = 0; l < elem.GetNP(); l++)
|
||||
{
|
||||
if (elem[l] == pi1) has1 = 1;
|
||||
if (elem[l] == pi2) has2 = 1;
|
||||
}
|
||||
|
||||
if (has1 && has2)
|
||||
{ // only once
|
||||
if (hasbothpoints.Contains (elnr))
|
||||
has1 = false;
|
||||
|
||||
if (has1)
|
||||
{
|
||||
hasbothpoints.Append (elnr);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (ElementIndex ei : hasbothpoints)
|
||||
{
|
||||
if (mesh[ei].GetType () != TET)
|
||||
return false;
|
||||
|
||||
if (mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
|
||||
return false;
|
||||
|
||||
|
||||
if ((mesh.ElementType(ei)) == FIXEDELEMENT)
|
||||
return false;
|
||||
|
||||
if(working_elements &&
|
||||
ei < working_elements->Size() &&
|
||||
!working_elements->Test(ei))
|
||||
return false;
|
||||
|
||||
if (mesh[ei].IsDeleted())
|
||||
return false;
|
||||
|
||||
if ((goal == OPT_LEGAL) &&
|
||||
mesh.LegalTet (mesh[ei]) &&
|
||||
CalcBad (mesh.Points(), mesh[ei], 0) < 1e3)
|
||||
return false;
|
||||
}
|
||||
|
||||
int nsuround = hasbothpoints.Size();
|
||||
int mattyp = mesh[hasbothpoints[0]].GetIndex();
|
||||
|
||||
if ( nsuround == 3 )
|
||||
{
|
||||
Element & elem = mesh[hasbothpoints[0]];
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elem[l] != pi1 && elem[l] != pi2)
|
||||
{
|
||||
pi4 = pi3;
|
||||
pi3 = elem[l];
|
||||
}
|
||||
|
||||
el31[0] = pi1;
|
||||
el31[1] = pi2;
|
||||
el31[2] = pi3;
|
||||
el31[3] = pi4;
|
||||
el31.SetIndex (mattyp);
|
||||
|
||||
if (WrongOrientation (mesh.Points(), el31))
|
||||
{
|
||||
Swap (pi3, pi4);
|
||||
el31[2] = pi3;
|
||||
el31[3] = pi4;
|
||||
}
|
||||
|
||||
pi5.Invalidate();
|
||||
for (int k = 0; k < 3; k++) // JS, 201212
|
||||
{
|
||||
const Element & elemk = mesh[hasbothpoints[k]];
|
||||
bool has1 = false;
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elemk[l] == pi4)
|
||||
has1 = true;
|
||||
if (has1)
|
||||
{
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elemk[l] != pi1 && elemk[l] != pi2 && elemk[l] != pi4)
|
||||
pi5 = elemk[l];
|
||||
}
|
||||
}
|
||||
|
||||
if (!pi5.IsValid())
|
||||
throw NgException("Illegal state observed in SwapImprove");
|
||||
|
||||
|
||||
el32[0] = pi1;
|
||||
el32[1] = pi2;
|
||||
el32[2] = pi4;
|
||||
el32[3] = pi5;
|
||||
el32.SetIndex (mattyp);
|
||||
|
||||
el33[0] = pi1;
|
||||
el33[1] = pi2;
|
||||
el33[2] = pi5;
|
||||
el33[3] = pi3;
|
||||
el33.SetIndex (mattyp);
|
||||
|
||||
bad1 = CalcBad (mesh.Points(), el31, 0) +
|
||||
CalcBad (mesh.Points(), el32, 0) +
|
||||
CalcBad (mesh.Points(), el33, 0);
|
||||
|
||||
el31.flags.illegal_valid = 0;
|
||||
el32.flags.illegal_valid = 0;
|
||||
el33.flags.illegal_valid = 0;
|
||||
|
||||
if (!mesh.LegalTet(el31) ||
|
||||
!mesh.LegalTet(el32) ||
|
||||
!mesh.LegalTet(el33))
|
||||
bad1 += 1e4;
|
||||
|
||||
el21[0] = pi3;
|
||||
el21[1] = pi4;
|
||||
el21[2] = pi5;
|
||||
el21[3] = pi2;
|
||||
el21.SetIndex (mattyp);
|
||||
|
||||
el22[0] = pi5;
|
||||
el22[1] = pi4;
|
||||
el22[2] = pi3;
|
||||
el22[3] = pi1;
|
||||
el22.SetIndex (mattyp);
|
||||
|
||||
bad2 = CalcBad (mesh.Points(), el21, 0) +
|
||||
CalcBad (mesh.Points(), el22, 0);
|
||||
|
||||
el21.flags.illegal_valid = 0;
|
||||
el22.flags.illegal_valid = 0;
|
||||
|
||||
if (!mesh.LegalTet(el21) ||
|
||||
!mesh.LegalTet(el22))
|
||||
bad2 += 1e4;
|
||||
|
||||
|
||||
if (goal == OPT_CONFORM && bad2 < 1e4)
|
||||
{
|
||||
INDEX_3 face(pi3, pi4, pi5);
|
||||
face.Sort();
|
||||
if (faces.Used(face))
|
||||
{
|
||||
// (*testout) << "3->2 swap, could improve conformity, bad1 = " << bad1
|
||||
// << ", bad2 = " << bad2 << endl;
|
||||
if (bad2 < 1e4)
|
||||
bad1 = 2 * bad2;
|
||||
}
|
||||
}
|
||||
|
||||
if (bad2 < bad1)
|
||||
{
|
||||
// (*mycout) << "3->2 " << flush;
|
||||
// (*testout) << "3->2 conversion" << endl;
|
||||
do_swap = true;
|
||||
if(check_only) return do_swap;
|
||||
|
||||
|
||||
/*
|
||||
(*testout) << "3->2 swap, old els = " << endl
|
||||
<< mesh[hasbothpoints[0]] << endl
|
||||
<< mesh[hasbothpoints[1]] << endl
|
||||
<< mesh[hasbothpoints[2]] << endl
|
||||
<< "new els = " << endl
|
||||
<< el21 << endl
|
||||
<< el22 << endl;
|
||||
*/
|
||||
|
||||
el21.flags.illegal_valid = 0;
|
||||
el22.flags.illegal_valid = 0;
|
||||
mesh[hasbothpoints[0]] = el21;
|
||||
mesh[hasbothpoints[1]] = el22;
|
||||
for (int l = 0; l < 4; l++)
|
||||
mesh[hasbothpoints[2]][l].Invalidate();
|
||||
mesh[hasbothpoints[2]].Delete();
|
||||
|
||||
elementsonnode.Add (pi4, hasbothpoints[1]);
|
||||
elementsonnode.Add (pi3, hasbothpoints[2]);
|
||||
|
||||
for (int k = 0; k < 2; k++)
|
||||
for (int l = 0; l < 4; l++)
|
||||
elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]);
|
||||
}
|
||||
}
|
||||
|
||||
if (nsuround == 4)
|
||||
{
|
||||
const Element & elem1 = mesh[hasbothpoints[0]];
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elem1[l] != pi1 && elem1[l] != pi2)
|
||||
{
|
||||
pi4 = pi3;
|
||||
pi3 = elem1[l];
|
||||
}
|
||||
|
||||
el1[0] = pi1; el1[1] = pi2;
|
||||
el1[2] = pi3; el1[3] = pi4;
|
||||
el1.SetIndex (mattyp);
|
||||
|
||||
if (WrongOrientation (mesh.Points(), el1))
|
||||
{
|
||||
Swap (pi3, pi4);
|
||||
el1[2] = pi3;
|
||||
el1[3] = pi4;
|
||||
}
|
||||
|
||||
pi5.Invalidate();
|
||||
for (int k = 0; k < 4; k++)
|
||||
{
|
||||
const Element & elem = mesh[hasbothpoints[k]];
|
||||
bool has1 = elem.PNums().Contains(pi4);
|
||||
if (has1)
|
||||
{
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elem[l] != pi1 && elem[l] != pi2 && elem[l] != pi4)
|
||||
pi5 = elem[l];
|
||||
}
|
||||
}
|
||||
|
||||
pi6.Invalidate();
|
||||
for (int k = 0; k < 4; k++)
|
||||
{
|
||||
const Element & elem = mesh[hasbothpoints[k]];
|
||||
bool has1 = elem.PNums().Contains(pi3);
|
||||
if (has1)
|
||||
{
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elem[l] != pi1 && elem[l] != pi2 && elem[l] != pi3)
|
||||
pi6 = elem[l];
|
||||
}
|
||||
}
|
||||
|
||||
el1[0] = pi1; el1[1] = pi2;
|
||||
el1[2] = pi3; el1[3] = pi4;
|
||||
el1.SetIndex (mattyp);
|
||||
|
||||
el2[0] = pi1; el2[1] = pi2;
|
||||
el2[2] = pi4; el2[3] = pi5;
|
||||
el2.SetIndex (mattyp);
|
||||
|
||||
el3[0] = pi1; el3[1] = pi2;
|
||||
el3[2] = pi5; el3[3] = pi6;
|
||||
el3.SetIndex (mattyp);
|
||||
|
||||
el4[0] = pi1; el4[1] = pi2;
|
||||
el4[2] = pi6; el4[3] = pi3;
|
||||
el4.SetIndex (mattyp);
|
||||
|
||||
bad1 = CalcBad (mesh.Points(), el1, 0) +
|
||||
CalcBad (mesh.Points(), el2, 0) +
|
||||
CalcBad (mesh.Points(), el3, 0) +
|
||||
CalcBad (mesh.Points(), el4, 0);
|
||||
|
||||
|
||||
el1.flags.illegal_valid = 0;
|
||||
el2.flags.illegal_valid = 0;
|
||||
el3.flags.illegal_valid = 0;
|
||||
el4.flags.illegal_valid = 0;
|
||||
|
||||
|
||||
if (goal != OPT_CONFORM)
|
||||
{
|
||||
if (!mesh.LegalTet(el1) ||
|
||||
!mesh.LegalTet(el2) ||
|
||||
!mesh.LegalTet(el3) ||
|
||||
!mesh.LegalTet(el4))
|
||||
bad1 += 1e4;
|
||||
}
|
||||
|
||||
el1[0] = pi3; el1[1] = pi5;
|
||||
el1[2] = pi2; el1[3] = pi4;
|
||||
el1.SetIndex (mattyp);
|
||||
|
||||
el2[0] = pi3; el2[1] = pi5;
|
||||
el2[2] = pi4; el2[3] = pi1;
|
||||
el2.SetIndex (mattyp);
|
||||
|
||||
el3[0] = pi3; el3[1] = pi5;
|
||||
el3[2] = pi1; el3[3] = pi6;
|
||||
el3.SetIndex (mattyp);
|
||||
|
||||
el4[0] = pi3; el4[1] = pi5;
|
||||
el4[2] = pi6; el4[3] = pi2;
|
||||
el4.SetIndex (mattyp);
|
||||
|
||||
bad2 = CalcBad (mesh.Points(), el1, 0) +
|
||||
CalcBad (mesh.Points(), el2, 0) +
|
||||
CalcBad (mesh.Points(), el3, 0) +
|
||||
CalcBad (mesh.Points(), el4, 0);
|
||||
|
||||
el1.flags.illegal_valid = 0;
|
||||
el2.flags.illegal_valid = 0;
|
||||
el3.flags.illegal_valid = 0;
|
||||
el4.flags.illegal_valid = 0;
|
||||
|
||||
if (goal != OPT_CONFORM)
|
||||
{
|
||||
if (!mesh.LegalTet(el1) ||
|
||||
!mesh.LegalTet(el2) ||
|
||||
!mesh.LegalTet(el3) ||
|
||||
!mesh.LegalTet(el4))
|
||||
bad2 += 1e4;
|
||||
}
|
||||
|
||||
|
||||
el1b[0] = pi4; el1b[1] = pi6;
|
||||
el1b[2] = pi3; el1b[3] = pi2;
|
||||
el1b.SetIndex (mattyp);
|
||||
|
||||
el2b[0] = pi4; el2b[1] = pi6;
|
||||
el2b[2] = pi2; el2b[3] = pi5;
|
||||
el2b.SetIndex (mattyp);
|
||||
|
||||
el3b[0] = pi4; el3b[1] = pi6;
|
||||
el3b[2] = pi5; el3b[3] = pi1;
|
||||
el3b.SetIndex (mattyp);
|
||||
|
||||
el4b[0] = pi4; el4b[1] = pi6;
|
||||
el4b[2] = pi1; el4b[3] = pi3;
|
||||
el4b.SetIndex (mattyp);
|
||||
|
||||
bad3 = CalcBad (mesh.Points(), el1b, 0) +
|
||||
CalcBad (mesh.Points(), el2b, 0) +
|
||||
CalcBad (mesh.Points(), el3b, 0) +
|
||||
CalcBad (mesh.Points(), el4b, 0);
|
||||
|
||||
el1b.flags.illegal_valid = 0;
|
||||
el2b.flags.illegal_valid = 0;
|
||||
el3b.flags.illegal_valid = 0;
|
||||
el4b.flags.illegal_valid = 0;
|
||||
|
||||
if (goal != OPT_CONFORM)
|
||||
{
|
||||
if (!mesh.LegalTet(el1b) ||
|
||||
!mesh.LegalTet(el2b) ||
|
||||
!mesh.LegalTet(el3b) ||
|
||||
!mesh.LegalTet(el4b))
|
||||
bad3 += 1e4;
|
||||
}
|
||||
|
||||
bool swap2, swap3;
|
||||
|
||||
if (goal != OPT_CONFORM)
|
||||
{
|
||||
swap2 = (bad2 < bad1) && (bad2 < bad3);
|
||||
swap3 = !swap2 && (bad3 < bad1);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (mesh.BoundaryEdge (pi3, pi5)) bad2 /= 1e6;
|
||||
if (mesh.BoundaryEdge (pi4, pi6)) bad3 /= 1e6;
|
||||
|
||||
swap2 = (bad2 < bad1) && (bad2 < bad3);
|
||||
swap3 = !swap2 && (bad3 < bad1);
|
||||
}
|
||||
|
||||
|
||||
if (swap2 || swap3)
|
||||
{
|
||||
// (*mycout) << "4->4 " << flush;
|
||||
do_swap = true;
|
||||
if(check_only) return do_swap;
|
||||
// (*testout) << "4->4 conversion" << "\n";
|
||||
/*
|
||||
(*testout) << "bad1 = " << bad1
|
||||
<< " bad2 = " << bad2
|
||||
<< " bad3 = " << bad3 << "\n";
|
||||
|
||||
(*testout) << "Points: " << pi1 << " " << pi2 << " " << pi3
|
||||
<< " " << pi4 << " " << pi5 << " " << pi6 << "\n";
|
||||
(*testout) << "Elements: "
|
||||
<< hasbothpoints.Get(1) << " "
|
||||
<< hasbothpoints.Get(2) << " "
|
||||
<< hasbothpoints.Get(3) << " "
|
||||
<< hasbothpoints.Get(4) << " " << "\n";
|
||||
*/
|
||||
|
||||
/*
|
||||
{
|
||||
int i1, j1;
|
||||
for (i1 = 1; i1 <= 4; i1++)
|
||||
{
|
||||
for (j1 = 1; j1 <= 4; j1++)
|
||||
(*testout) << volelements.Get(hasbothpoints.Get(i1)).PNum(j1)
|
||||
<< " ";
|
||||
(*testout) << "\n";
|
||||
}
|
||||
}
|
||||
*/
|
||||
}
|
||||
|
||||
|
||||
if (swap2)
|
||||
{
|
||||
// (*mycout) << "bad1 = " << bad1 << " bad2 = " << bad2 << "\n";
|
||||
|
||||
|
||||
/*
|
||||
(*testout) << "4->4 swap A, old els = " << endl
|
||||
<< mesh[hasbothpoints[0]] << endl
|
||||
<< mesh[hasbothpoints[1]] << endl
|
||||
<< mesh[hasbothpoints[2]] << endl
|
||||
<< mesh[hasbothpoints[3]] << endl
|
||||
<< "new els = " << endl
|
||||
<< el1 << endl
|
||||
<< el2 << endl
|
||||
<< el3 << endl
|
||||
<< el4 << endl;
|
||||
*/
|
||||
|
||||
|
||||
|
||||
el1.flags.illegal_valid = 0;
|
||||
el2.flags.illegal_valid = 0;
|
||||
el3.flags.illegal_valid = 0;
|
||||
el4.flags.illegal_valid = 0;
|
||||
|
||||
mesh[hasbothpoints[0]] = el1;
|
||||
mesh[hasbothpoints[1]] = el2;
|
||||
mesh[hasbothpoints[2]] = el3;
|
||||
mesh[hasbothpoints[3]] = el4;
|
||||
|
||||
for (int k = 0; k < 4; k++)
|
||||
for (int l = 0; l < 4; l++)
|
||||
elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]);
|
||||
}
|
||||
else if (swap3)
|
||||
{
|
||||
// (*mycout) << "bad1 = " << bad1 << " bad3 = " << bad3 << "\n";
|
||||
el1b.flags.illegal_valid = 0;
|
||||
el2b.flags.illegal_valid = 0;
|
||||
el3b.flags.illegal_valid = 0;
|
||||
el4b.flags.illegal_valid = 0;
|
||||
|
||||
|
||||
/*
|
||||
(*testout) << "4->4 swap A, old els = " << endl
|
||||
<< mesh[hasbothpoints[0]] << endl
|
||||
<< mesh[hasbothpoints[1]] << endl
|
||||
<< mesh[hasbothpoints[2]] << endl
|
||||
<< mesh[hasbothpoints[3]] << endl
|
||||
<< "new els = " << endl
|
||||
<< el1b << endl
|
||||
<< el2b << endl
|
||||
<< el3b << endl
|
||||
<< el4b << endl;
|
||||
*/
|
||||
|
||||
|
||||
mesh[hasbothpoints[0]] = el1b;
|
||||
mesh[hasbothpoints[1]] = el2b;
|
||||
mesh[hasbothpoints[2]] = el3b;
|
||||
mesh[hasbothpoints[3]] = el4b;
|
||||
|
||||
for (int k = 0; k < 4; k++)
|
||||
for (int l = 0; l < 4; l++)
|
||||
elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]);
|
||||
}
|
||||
}
|
||||
|
||||
// if (goal == OPT_QUALITY)
|
||||
if (nsuround >= 5)
|
||||
{
|
||||
Element hel(TET);
|
||||
|
||||
NgArrayMem<PointIndex, 50> suroundpts(nsuround);
|
||||
NgArrayMem<bool, 50> tetused(nsuround);
|
||||
|
||||
Element & elem = mesh[hasbothpoints[0]];
|
||||
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elem[l] != pi1 && elem[l] != pi2)
|
||||
{
|
||||
pi4 = pi3;
|
||||
pi3 = elem[l];
|
||||
}
|
||||
|
||||
hel[0] = pi1;
|
||||
hel[1] = pi2;
|
||||
hel[2] = pi3;
|
||||
hel[3] = pi4;
|
||||
hel.SetIndex (mattyp);
|
||||
|
||||
if (WrongOrientation (mesh.Points(), hel))
|
||||
{
|
||||
Swap (pi3, pi4);
|
||||
hel[2] = pi3;
|
||||
hel[3] = pi4;
|
||||
}
|
||||
|
||||
|
||||
// suroundpts.SetSize (nsuround);
|
||||
suroundpts = PointIndex::INVALID;
|
||||
suroundpts[0] = pi3;
|
||||
suroundpts[1] = pi4;
|
||||
|
||||
tetused = false;
|
||||
tetused[0] = true;
|
||||
|
||||
for (int l = 2; l < nsuround; l++)
|
||||
{
|
||||
PointIndex oldpi = suroundpts[l-1];
|
||||
PointIndex newpi;
|
||||
newpi.Invalidate();
|
||||
|
||||
for (int k = 0; k < nsuround && !newpi.IsValid(); k++)
|
||||
if (!tetused[k])
|
||||
{
|
||||
const Element & nel = mesh[hasbothpoints[k]];
|
||||
for (int k2 = 0; k2 < 4 && !newpi.IsValid(); k2++)
|
||||
if (nel[k2] == oldpi)
|
||||
{
|
||||
newpi =
|
||||
nel[0] + nel[1] + nel[2] + nel[3]
|
||||
- pi1 - pi2 - oldpi;
|
||||
|
||||
tetused[k] = true;
|
||||
suroundpts[l] = newpi;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
bad1 = 0;
|
||||
for (int k = 0; k < nsuround; k++)
|
||||
{
|
||||
hel[0] = pi1;
|
||||
hel[1] = pi2;
|
||||
hel[2] = suroundpts[k];
|
||||
hel[3] = suroundpts[(k+1) % nsuround];
|
||||
hel.SetIndex (mattyp);
|
||||
|
||||
bad1 += CalcBad (mesh.Points(), hel, 0);
|
||||
}
|
||||
|
||||
// (*testout) << "nsuround = " << nsuround << " bad1 = " << bad1 << endl;
|
||||
|
||||
|
||||
int bestl = -1;
|
||||
int confface = -1;
|
||||
int confedge = -1;
|
||||
double badopt = bad1;
|
||||
|
||||
for (int l = 0; l < nsuround; l++)
|
||||
{
|
||||
bad2 = 0;
|
||||
|
||||
for (int k = l+1; k <= nsuround + l - 2; k++)
|
||||
{
|
||||
hel[0] = suroundpts[l];
|
||||
hel[1] = suroundpts[k % nsuround];
|
||||
hel[2] = suroundpts[(k+1) % nsuround];
|
||||
hel[3] = pi2;
|
||||
|
||||
bad2 += CalcBad (mesh.Points(), hel, 0);
|
||||
hel.flags.illegal_valid = 0;
|
||||
if (!mesh.LegalTet(hel)) bad2 += 1e4;
|
||||
|
||||
hel[2] = suroundpts[k % nsuround];
|
||||
hel[1] = suroundpts[(k+1) % nsuround];
|
||||
hel[3] = pi1;
|
||||
|
||||
bad2 += CalcBad (mesh.Points(), hel, 0);
|
||||
|
||||
hel.flags.illegal_valid = 0;
|
||||
if (!mesh.LegalTet(hel)) bad2 += 1e4;
|
||||
}
|
||||
// (*testout) << "bad2," << l << " = " << bad2 << endl;
|
||||
|
||||
if ( bad2 < badopt )
|
||||
{
|
||||
bestl = l;
|
||||
badopt = bad2;
|
||||
}
|
||||
|
||||
|
||||
if (goal == OPT_CONFORM)
|
||||
// (bad2 <= 100 * bad1 || bad2 <= 1e6))
|
||||
{
|
||||
bool nottoobad =
|
||||
(bad2 <= bad1) ||
|
||||
(bad2 <= 100 * bad1 && bad2 <= 1e18) ||
|
||||
(bad2 <= 1e8);
|
||||
|
||||
for (int k = l+1; k <= nsuround + l - 2; k++)
|
||||
{
|
||||
INDEX_3 hi3(suroundpts[l],
|
||||
suroundpts[k % nsuround],
|
||||
suroundpts[(k+1) % nsuround]);
|
||||
hi3.Sort();
|
||||
if (faces.Used(hi3))
|
||||
{
|
||||
// (*testout) << "could improve face conformity, bad1 = " << bad1
|
||||
// << ", bad 2 = " << bad2 << ", nottoobad = " << nottoobad << endl;
|
||||
if (nottoobad)
|
||||
confface = l;
|
||||
}
|
||||
}
|
||||
|
||||
for (int k = l+2; k <= nsuround+l-2; k++)
|
||||
{
|
||||
if (mesh.BoundaryEdge (suroundpts[l],
|
||||
suroundpts[k % nsuround]))
|
||||
{
|
||||
/*
|
||||
*testout << "could improve edge conformity, bad1 = " << bad1
|
||||
<< ", bad 2 = " << bad2 << ", nottoobad = " << nottoobad << endl;
|
||||
*/
|
||||
if (nottoobad)
|
||||
confedge = l;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (confedge != -1)
|
||||
bestl = confedge;
|
||||
if (confface != -1)
|
||||
bestl = confface;
|
||||
|
||||
if (bestl != -1)
|
||||
{
|
||||
// (*mycout) << nsuround << "->" << 2 * (nsuround-2) << " " << flush;
|
||||
do_swap = true;
|
||||
if(check_only) return do_swap;
|
||||
|
||||
for (int k = bestl+1; k <= nsuround + bestl - 2; k++)
|
||||
{
|
||||
int k1;
|
||||
|
||||
hel[0] = suroundpts[bestl];
|
||||
hel[1] = suroundpts[k % nsuround];
|
||||
hel[2] = suroundpts[(k+1) % nsuround];
|
||||
hel[3] = pi2;
|
||||
hel.flags.illegal_valid = 0;
|
||||
|
||||
/*
|
||||
(*testout) << nsuround << "-swap, new el,top = "
|
||||
<< hel << endl;
|
||||
*/
|
||||
mesh.AddVolumeElement (hel);
|
||||
|
||||
for (k1 = 0; k1 < 4; k1++)
|
||||
elementsonnode.Add (hel[k1], mesh.GetNE()-1);
|
||||
|
||||
|
||||
hel[2] = suroundpts[k % nsuround];
|
||||
hel[1] = suroundpts[(k+1) % nsuround];
|
||||
hel[3] = pi1;
|
||||
|
||||
/*
|
||||
(*testout) << nsuround << "-swap, new el,bot = "
|
||||
<< hel << endl;
|
||||
*/
|
||||
|
||||
mesh.AddVolumeElement (hel);
|
||||
|
||||
for (k1 = 0; k1 < 4; k1++)
|
||||
elementsonnode.Add (hel[k1], mesh.GetNE()-1);
|
||||
}
|
||||
|
||||
for (int k = 0; k < nsuround; k++)
|
||||
{
|
||||
Element & rel = mesh[hasbothpoints[k]];
|
||||
/*
|
||||
(*testout) << nsuround << "-swap, old el = "
|
||||
<< rel << endl;
|
||||
*/
|
||||
rel.Delete();
|
||||
for (int k1 = 0; k1 < 4; k1++)
|
||||
rel[k1].Invalidate();
|
||||
}
|
||||
}
|
||||
}
|
||||
return do_swap;
|
||||
}
|
||||
|
||||
void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
const NgBitArray * working_elements)
|
||||
{
|
||||
static Timer t("MeshOptimize3d::SwapImprove"); RegionTimer reg(t);
|
||||
static Timer tloop("MeshOptimize3d::SwapImprove loop");
|
||||
|
||||
// return SwapImproveSequential(mesh, goal, working_elements);
|
||||
|
||||
int cnt = 0;
|
||||
|
||||
int np = mesh.GetNP();
|
||||
int ne = mesh.GetNE();
|
||||
|
||||
mesh.BoundaryEdge (1,2); // ensure the boundary-elements table is built
|
||||
|
||||
// contains at least all elements at node
|
||||
TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np);
|
||||
|
||||
NgArray<ElementIndex> hasbothpoints;
|
||||
|
||||
PrintMessage (3, "SwapImprove ");
|
||||
(*testout) << "\n" << "Start SwapImprove" << endl;
|
||||
|
||||
const char * savetask = multithread.task;
|
||||
multithread.task = "Swap Improve";
|
||||
|
||||
INDEX_3_HASHTABLE<int> faces(mesh.GetNOpenElements()/3 + 2);
|
||||
if (goal == OPT_CONFORM)
|
||||
{
|
||||
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
|
||||
{
|
||||
const Element2d & hel = mesh.OpenElement(i);
|
||||
INDEX_3 face(hel[0], hel[1], hel[2]);
|
||||
face.Sort();
|
||||
faces.Set (face, 1);
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate total badness
|
||||
if (goal == OPT_QUALITY)
|
||||
{
|
||||
double bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
|
||||
(*testout) << "Total badness = " << bad1 << endl;
|
||||
}
|
||||
|
||||
// find elements on node
|
||||
for (ElementIndex ei = 0; ei < ne; ei++)
|
||||
for (PointIndex pi : mesh[ei].PNums())
|
||||
elementsonnode.Add (pi, ei);
|
||||
|
||||
Array<std::tuple<PointIndex,PointIndex>> edges;
|
||||
BuildEdgeList(mesh, elementsonnode, edges);
|
||||
|
||||
Array<int> candidate_edges(edges.Size());
|
||||
std::atomic<int> improvement_counter(0);
|
||||
|
||||
tloop.Start();
|
||||
|
||||
ParallelForRange(Range(edges), [&] (auto myrange)
|
||||
{
|
||||
for(auto i : myrange)
|
||||
{
|
||||
if (multithread.terminate)
|
||||
break;
|
||||
|
||||
auto [pi0, pi1] = edges[i];
|
||||
if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, true))
|
||||
candidate_edges[improvement_counter++] = i;
|
||||
}
|
||||
});
|
||||
|
||||
auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
|
||||
QuickSort(edges_with_improvement);
|
||||
|
||||
for(auto ei : edges_with_improvement)
|
||||
{
|
||||
auto [pi0,pi1] = edges[ei];
|
||||
if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, false))
|
||||
cnt++;
|
||||
}
|
||||
|
||||
tloop.Stop();
|
||||
|
||||
PrintMessage (5, cnt, " swaps performed");
|
||||
|
||||
mesh.Compress ();
|
||||
|
||||
multithread.task = savetask;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
@ -26,8 +26,12 @@ public:
|
||||
void CombineImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
|
||||
|
||||
void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
|
||||
|
||||
bool SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, const NgBitArray * working_elements, TABLE<ElementIndex,PointIndex::BASE> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false);
|
||||
void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
|
||||
const NgBitArray * working_elements = NULL);
|
||||
void SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
|
||||
const NgBitArray * working_elements = NULL);
|
||||
void SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
|
||||
const NgBitArray * working_elements = NULL,
|
||||
const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL);
|
||||
|
||||
@ -506,7 +506,7 @@ namespace netgen
|
||||
cout << "set debugflag" << endl;
|
||||
}
|
||||
|
||||
if (debugparam.haltlargequalclass && qualclass > 50)
|
||||
if (debugparam.haltlargequalclass && qualclass == 50)
|
||||
debugflag = 1;
|
||||
|
||||
// problem recognition !
|
||||
|
||||
@ -19,6 +19,9 @@ int chartdebug = 0;
|
||||
|
||||
void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, const STLParameters& stlparam)
|
||||
{
|
||||
static Timer t("makeatlas"); RegionTimer reg(t);
|
||||
static Timer tinner("find innner chart");
|
||||
static Timer touter("find outer chart");
|
||||
// int timer1 = NgProfiler::CreateTimer ("makeatlas");
|
||||
/*
|
||||
int timerb = NgProfiler::CreateTimer ("makeatlas - begin");
|
||||
@ -183,7 +186,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
|
||||
// NgProfiler::StopTimer (timerb);
|
||||
// NgProfiler::StartTimer (timer2);
|
||||
|
||||
|
||||
tinner.Start();
|
||||
while (changed)
|
||||
{
|
||||
changed = false;
|
||||
@ -197,7 +200,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
|
||||
{
|
||||
for (int j = 1; j <= NONeighbourTrigs(i); j++)
|
||||
{
|
||||
int nt = NeighbourTrig(i,j);
|
||||
STLTrigId nt = NeighbourTrig(i,j);
|
||||
// *testout << "check trig " << nt << endl;
|
||||
STLPointId np1, np2;
|
||||
GetTriangle(i).GetNeighbourPoints(GetTriangle(nt),np1,np2);
|
||||
@ -316,7 +319,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
tinner.Stop();
|
||||
|
||||
innerchartpts.SetSize(innerchartpoints.Size());
|
||||
for (size_t i = 0; i < innerchartpoints.Size(); i++)
|
||||
innerchartpts[i] = GetPoint(innerchartpoints[i]);
|
||||
@ -329,8 +333,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
|
||||
|
||||
// chartbound.Clear();
|
||||
// warum, ic-bound auf edge macht Probleme js ???
|
||||
|
||||
|
||||
|
||||
touter.Start();
|
||||
outermark[starttrig] = chartnum;
|
||||
//chart->AddOuterTrig(starttrig);
|
||||
changed = true;
|
||||
@ -500,7 +504,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
touter.Stop();
|
||||
// NgProfiler::StopTimer (timer3);
|
||||
// NgProfiler::StartTimer (timere);
|
||||
// NgProfiler::StartTimer (timere1);
|
||||
|
||||
@ -867,6 +867,7 @@ void STLBoundary :: AddOrDelSegment(const STLBoundarySeg & seg)
|
||||
|
||||
void STLBoundary ::AddTriangle(const STLTriangle & t)
|
||||
{
|
||||
// static Timer timer("STLBoundary::AddTriangle"); RegionTimer reg(timer);
|
||||
// static int timer_old = NgProfiler::CreateTimer ("STLChart::AddTriangle_old");
|
||||
// static int timer_new = NgProfiler::CreateTimer ("STLChart::AddTriangle_new");
|
||||
|
||||
@ -1009,9 +1010,17 @@ void STLBoundary ::AddTriangle(const STLTriangle & t)
|
||||
|
||||
INDEX_2 op(seg[1], seg[0]);
|
||||
if (boundary_ht.Used(op))
|
||||
boundary_ht.Delete(op);
|
||||
{
|
||||
boundary_ht.Delete(op);
|
||||
if (searchtree)
|
||||
searchtree->DeleteElement(op);
|
||||
}
|
||||
else
|
||||
boundary_ht[seg] = bseg;
|
||||
{
|
||||
boundary_ht[seg] = bseg;
|
||||
if (searchtree)
|
||||
searchtree->Insert (bseg.BoundingBox(), seg);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -1211,54 +1220,32 @@ bool STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3
|
||||
|
||||
void STLBoundary :: BuildSearchTree()
|
||||
{
|
||||
// static int timer = NgProfiler::CreateTimer ("BuildSearchTree");
|
||||
// NgProfiler::RegionTimer reg(timer);
|
||||
|
||||
delete searchtree;
|
||||
|
||||
/*
|
||||
Box<2> box2d(Box<2>::EMPTY_BOX);
|
||||
|
||||
int nseg = NOSegments();
|
||||
for (int j = 1; j <= nseg; j++)
|
||||
{
|
||||
const STLBoundarySeg & seg = GetSegment(j);
|
||||
if (seg.IsSmoothEdge()) continue;
|
||||
box2d.Add(seg.BoundingBox().PMin());
|
||||
box2d.Add(seg.BoundingBox().PMax());
|
||||
}
|
||||
|
||||
searchtree = new BoxTree<2> (box2d);
|
||||
|
||||
for (int j = 1; j <= nseg; j++)
|
||||
{
|
||||
const STLBoundarySeg & seg = GetSegment(j);
|
||||
if (seg.IsSmoothEdge()) continue;
|
||||
searchtree -> Insert (seg.BoundingBox(), j);
|
||||
}
|
||||
*/
|
||||
|
||||
Box<2> box2d(Box<2>::EMPTY_BOX);
|
||||
Box<3> box3d = geometry->GetBoundingBox();
|
||||
for (size_t i = 0; i < 8; i++)
|
||||
box2d.Add ( chart->Project2d (box3d.GetPointNr(i)));
|
||||
searchtree = new BoxTree<2,INDEX_2> (box2d);
|
||||
|
||||
// comment to enable searchtree:
|
||||
// searchtree = new BoxTree<2,INDEX_2> (box2d);
|
||||
searchtree = nullptr;
|
||||
}
|
||||
|
||||
void STLBoundary :: DeleteSearchTree()
|
||||
{
|
||||
// static int timer = NgProfiler::CreateTimer ("DeleteSearchTree");
|
||||
// NgProfiler::RegionTimer reg(timer);
|
||||
|
||||
delete searchtree;
|
||||
searchtree = nullptr;
|
||||
}
|
||||
|
||||
|
||||
// checks, whether 2d projection intersects
|
||||
bool STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
|
||||
const Vec3d& sn)
|
||||
{
|
||||
// static int timerquick = NgProfiler::CreateTimer ("TestSegChartNV-searchtree");
|
||||
// static int timer = NgProfiler::CreateTimer ("TestSegChartNV");
|
||||
// static int timerquick = NgProfiler::CreateTimer ("TestSegChartNV-searchtree");
|
||||
// static Timer timer("TestSegChartNV"); RegionTimer reg(timer);
|
||||
|
||||
|
||||
Point<2> p2d1 = chart->Project2d (p1);
|
||||
Point<2> p2d2 = chart->Project2d (p2);
|
||||
@ -1272,26 +1259,12 @@ bool STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
|
||||
double eps = 1e-6;
|
||||
bool ok = true;
|
||||
|
||||
/*
|
||||
static long int cnt = 0;
|
||||
static long int totnseg = 0;
|
||||
totnseg += nseg;
|
||||
cnt++;
|
||||
if ( (cnt % 100000) == 0)
|
||||
cout << "avg nseg = " << double(totnseg)/cnt << endl;
|
||||
*/
|
||||
|
||||
// TODO: fix searchtree update
|
||||
if (false)
|
||||
if (searchtree)
|
||||
{
|
||||
// NgProfiler::RegionTimer reg(timerquick);
|
||||
|
||||
NgArrayMem<INDEX_2,100> pis;
|
||||
searchtree -> GetIntersecting (box2d.PMin(), box2d.PMax(), pis);
|
||||
|
||||
for (auto i2 : pis)
|
||||
{
|
||||
// const STLBoundarySeg & seg = GetSegment(j);
|
||||
const STLBoundarySeg & seg = boundary_ht[i2];
|
||||
|
||||
if (seg.IsSmoothEdge()) continue;
|
||||
@ -1311,7 +1284,6 @@ bool STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
|
||||
|
||||
if(!err && ((on1 && in2) || (on2 && in1)))
|
||||
{
|
||||
|
||||
ok = false;
|
||||
break;
|
||||
}
|
||||
@ -1320,7 +1292,6 @@ bool STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
|
||||
|
||||
else
|
||||
{
|
||||
// NgProfiler::RegionTimer reg(timer);
|
||||
for(auto [i2, seg] : boundary_ht)
|
||||
{
|
||||
if (seg.IsSmoothEdge()) continue;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user