ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-02-11 21:03:07 +05:00
Code Issues Packages Projects Releases Wiki Activity

PointIndex in improve2gen

Browse Source
This commit is contained in:
Joachim Schoeberl 2024-12-28 10:35:18 +01:00
parent 2291221719
commit 10a56a9e86
4 changed files with 133 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
libsrc/core/array.hpp
View File

@ -554,6 +554,13 @@ namespace ngcore
// { return CArray<T> (data+pos); }
NETGEN_INLINE T * operator+ (size_t pos) const { return data+pos; }
/// access first element. check by macro NETGEN_CHECK_RANGE
T & First () const
{
NETGEN_CHECK_RANGE(0,0,size);
return data[0];
}
/// access last element. check by macro NETGEN_CHECK_RANGE
T & Last () const
{

213
libsrc/meshing/improve2gen.cpp
View File

@ -11,9 +11,9 @@ namespace netgen
public:
NgArray<Element2d> oldels;
NgArray<Element2d> newels;
NgArray<INDEX_2> deledges;
NgArray<int> incelsonnode;
NgArray<int> reused;
NgArray<PointIndices<2>> deledges;
Array<int,PointIndex> incelsonnode;
Array<int,PointIndex> reused;
int bonus;
int onp;
};
@ -55,128 +55,137 @@ namespace netgen
NgArray<ImprovementRule*> rules;
NgArray<SurfaceElementIndex> elmap;
NgArray<int> elrot;
NgArray<PointIndex> pmap;
NgArray<PointGeomInfo> pgi;
Array<PointIndex,PointIndex> pmap;
Array<PointGeomInfo,PointIndex> pgi;
int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
ImprovementRule * r1;
PointIndex p1 = IndexBASE<PointIndex>();
PointIndex p2 = p1+1;
PointIndex p3 = p2+1;
PointIndex p4 = p3+1;
PointIndex p5 = p4+1;
PointIndex p6 = p5+1;
PointIndex p7 = p6+1;
// 2 triangles to quad
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3));
r1->oldels.Append (Element2d (3, 2, 4));
r1->newels.Append (Element2d (1, 2, 4, 3));
r1->deledges.Append (INDEX_2 (2,3));
r1->oldels.Append (Element2d (p1, p2, p3));
r1->oldels.Append (Element2d (p3, p2, p4));
r1->newels.Append (Element2d (p1, p2, p4, p3));
r1->deledges.Append ( { p2, p3 } );
r1->onp = 4;
r1->bonus = 2;
rules.Append (r1);
// 2 quad to 1 quad
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (4, 3, 2, 5));
r1->newels.Append (Element2d (1, 2, 5, 4));
r1->deledges.Append (INDEX_2 (2, 3));
r1->deledges.Append (INDEX_2 (3, 4));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p4, p3, p2, p5));
r1->newels.Append (Element2d (p1, p2, p5, p4));
r1->deledges.Append ( { p2, p3 } );
r1->deledges.Append ( { p3, p4 } );
r1->onp = 5;
r1->bonus = 0;
rules.Append (r1);
// swap quads
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (3, 2, 5, 6));
r1->newels.Append (Element2d (1, 6, 3, 4));
r1->newels.Append (Element2d (1, 2, 5, 6));
r1->deledges.Append (INDEX_2 (2, 3));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p3, p2, p5, p6));
r1->newels.Append (Element2d (p1, p6, p3, p4));
r1->newels.Append (Element2d (p1, p2, p5, p6));
r1->deledges.Append ( { p2, p3 } );
r1->onp = 6;
r1->bonus = 0;
rules.Append (r1);
// three quads to 2
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (2, 5, 6, 3));
r1->oldels.Append (Element2d (3, 6, 7, 4));
r1->newels.Append (Element2d (1, 2, 5, 4));
r1->newels.Append (Element2d (4, 5, 6, 7));
r1->deledges.Append (INDEX_2 (2, 3));
r1->deledges.Append (INDEX_2 (3, 4));
r1->deledges.Append (INDEX_2 (3, 6));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p2, p5, p6, p3));
r1->oldels.Append (Element2d (p3, p6, p7, p4));
r1->newels.Append (Element2d (p1, p2, p5, p4));
r1->newels.Append (Element2d (p4, p5, p6, p7));
r1->deledges.Append ( { p2, p3 } );
r1->deledges.Append ( { p3, p4 } );
r1->deledges.Append ( { p3, p6 } );
r1->onp = 7;
r1->bonus = -1;
rules.Append (r1);
// quad + 2 connected trigs to quad
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (2, 5, 3));
r1->oldels.Append (Element2d (3, 5, 4));
r1->newels.Append (Element2d (1, 2, 5, 4));
r1->deledges.Append (INDEX_2 (2, 3));
r1->deledges.Append (INDEX_2 (3, 4));
r1->deledges.Append (INDEX_2 (3, 5));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p2, p5, p3));
r1->oldels.Append (Element2d (p3, p5, p4));
r1->newels.Append (Element2d (p1, p2, p5, p4));
r1->deledges.Append ( { p2, p3 } );
r1->deledges.Append ( { p3, p4 } );
r1->deledges.Append ( { p3, p5 } );
r1->onp = 5;
r1->bonus = 0;
rules.Append (r1);
// quad + 2 non-connected trigs to quad (a and b)
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (2, 6, 3));
r1->oldels.Append (Element2d (1, 4, 5));
r1->newels.Append (Element2d (1, 3, 4, 5));
r1->newels.Append (Element2d (1, 2, 6, 3));
r1->deledges.Append (INDEX_2 (1, 4));
r1->deledges.Append (INDEX_2 (2, 3));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p2, p6, p3));
r1->oldels.Append (Element2d (p1, p4, p5));
r1->newels.Append (Element2d (p1, p3, p4, p5));
r1->newels.Append (Element2d (p1, p2, p6, p3));
r1->deledges.Append ( { p1, p4 } );
r1->deledges.Append ( { p2, p3 } );
r1->onp = 6;
r1->bonus = 0;
rules.Append (r1);
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (2, 6, 3));
r1->oldels.Append (Element2d (1, 4, 5));
r1->newels.Append (Element2d (1, 2, 4, 5));
r1->newels.Append (Element2d (4, 2, 6, 3));
r1->deledges.Append (INDEX_2 (1, 4));
r1->deledges.Append (INDEX_2 (2, 3));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p2, p6, p3));
r1->oldels.Append (Element2d (p1, p4, p5));
r1->newels.Append (Element2d (p1, p2, p4, p5));
r1->newels.Append (Element2d (p4, p2, p6, p3));
r1->deledges.Append ( { p1, p4 } );
r1->deledges.Append ( { p2, p3 } );
r1->onp = 6;
r1->bonus = 0;
rules.Append (r1);
// two quad + trig -> one quad + trig
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (2, 5, 6, 3));
r1->oldels.Append (Element2d (4, 3, 6));
r1->newels.Append (Element2d (1, 2, 6, 4));
r1->newels.Append (Element2d (2, 5, 6));
r1->deledges.Append (INDEX_2 (2, 3));
r1->deledges.Append (INDEX_2 (3, 4));
r1->deledges.Append (INDEX_2 (3, 6));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p2, p5, p6, p3));
r1->oldels.Append (Element2d (p4, p3, p6));
r1->newels.Append (Element2d (p1, p2, p6, p4));
r1->newels.Append (Element2d (p2, p5, p6));
r1->deledges.Append ( { p2, p3 } );
r1->deledges.Append ( { p3, p4 } );
r1->deledges.Append ( { p3, p6 } );
r1->onp = 6;
r1->bonus = -1;
rules.Append (r1);
// swap quad + trig (a and b)
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (2, 5, 3));
r1->newels.Append (Element2d (2, 5, 3, 4));
r1->newels.Append (Element2d (1, 2, 4));
r1->deledges.Append (INDEX_2 (2, 3));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p2, p5, p3));
r1->newels.Append (Element2d (p2, p5, p3, p4));
r1->newels.Append (Element2d (p1, p2, p4));
r1->deledges.Append ( { p2, p3 } );
r1->onp = 5;
r1->bonus = 0;
rules.Append (r1);
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (2, 5, 3));
r1->newels.Append (Element2d (1, 2, 5, 3));
r1->newels.Append (Element2d (1, 3, 4));
r1->deledges.Append (INDEX_2 (2, 3));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p2, p5, p3));
r1->newels.Append (Element2d (p1, p2, p5, p3));
r1->newels.Append (Element2d (p1, p3, p4));
r1->deledges.Append ( { p2, p3 } );
r1->onp = 5;
r1->bonus = 0;
rules.Append (r1);
@ -184,12 +193,12 @@ namespace netgen
// 2 quads to quad + 2 trigs
r1 = new ImprovementRule;
r1->oldels.Append (Element2d (1, 2, 3, 4));
r1->oldels.Append (Element2d (3, 2, 5, 6));
r1->newels.Append (Element2d (1, 5, 6, 4));
r1->newels.Append (Element2d (1, 2, 5));
r1->newels.Append (Element2d (4, 6, 3));
r1->deledges.Append (INDEX_2 (2, 3));
r1->oldels.Append (Element2d (p1, p2, p3, p4));
r1->oldels.Append (Element2d (p3, p2, p5, p6));
r1->newels.Append (Element2d (p1, p5, p6, p4));
r1->newels.Append (Element2d (p1, p2, p5));
r1->newels.Append (Element2d (p4, p6, p3));
r1->deledges.Append ( { p2, p3 } );
r1->onp = 6;
r1->bonus = 0;
// rules.Append (r1);
@ -203,18 +212,21 @@ namespace netgen
mapped = 0;
for (int ri = 0; ri < rules.Size(); ri++)
{
ImprovementRule & rule = *rules[ri];
rule.incelsonnode.SetSize (rule.onp);
rule.reused.SetSize (rule.onp);
/*
for (int j = 0; j < rule.onp; j++)
{
rule.incelsonnode[j] = 0;
rule.reused[j] = 0;
}
*/
rule.incelsonnode = 0;
rule.reused = 0;
/*
for (int j = 1; j <= rule.oldels.Size(); j++)
@ -226,15 +238,15 @@ namespace netgen
*/
for (const Element2d & el : rule.oldels)
for (PointIndex pi : el.PNums())
rule.incelsonnode[pi-IndexBASE<PointIndex>()]--;
rule.incelsonnode[pi]--;
for (int j = 1; j <= rule.newels.Size(); j++)
{
const Element2d & el = rule.newels.Elem(j);
for (int k = 1; k <= el.GetNP(); k++)
{
rule.incelsonnode.Elem(el.PNum(k))++;
rule.reused.Elem(el.PNum(k)) = 1;
rule.incelsonnode[el.PNum(k)]++;
rule.reused[el.PNum(k)] = 1;
}
}
}
@ -242,8 +254,8 @@ namespace netgen
TABLE<int,PointIndex::BASE> elonnode(np);
NgArray<int,PointIndex::BASE> nelonnode(np);
DynamicTable<int,PointIndex> elonnode(np);
Array<int,PointIndex> nelonnode(np);
TABLE<SurfaceElementIndex> nbels(ne);
nelonnode = -4;
@ -314,13 +326,13 @@ namespace netgen
if (el0.IsDeleted()) continue;
if (el0.GetNP() != rel0.GetNP()) continue;
pmap = PointIndex (-1);
// pmap = PointIndex (-1);
for (auto & p : pmap) p.Invalidate();
for (int k = 0; k < el0.GetNP(); k++)
{
pmap.Elem(rel0[k]) = el0.PNumMod(k+elrot[0]+1);
pgi.Elem(rel0[k]) = el0.GeomInfoPiMod(k+elrot[0]+1);
pmap[rel0[k]] = el0.PNumMod(k+elrot[0]+1);
pgi[rel0[k]] = el0.GeomInfoPiMod(k+elrot[0]+1);
}
ok = 1;
@ -342,16 +354,16 @@ namespace netgen
possible = 1;
for (int k = 0; k < rel.GetNP(); k++)
if (pmap.Elem(rel[k]).IsValid() &&
pmap.Elem(rel[k]) != el.PNumMod (k+elrot[i]+1))
if (pmap[rel[k]].IsValid() &&
pmap[rel[k]] != el.PNumMod (k+elrot[i]+1))
possible = 0;
if (possible)
{
for (int k = 0; k < el.GetNP(); k++)
{
pmap.Elem(rel[k]) = el.PNumMod(k+elrot[i]+1);
pgi.Elem(rel[k]) = el.GeomInfoPiMod(k+elrot[i]+1);
pmap[rel[k]] = el.PNumMod(k+elrot[i]+1);
pgi[rel[k]] = el.GeomInfoPiMod(k+elrot[i]+1);
}
break;
}
@ -370,8 +382,8 @@ namespace netgen
for(int i=0; ok && i<rule.deledges.Size(); i++)
{
ok = !mesh.IsSegment(pmap.Elem(rule.deledges[i].I1()),
pmap.Elem(rule.deledges[i].I2()));
ok = !mesh.IsSegment(pmap[rule.deledges[i][0]],
pmap[rule.deledges[i][1]]);
}
@ -382,15 +394,15 @@ namespace netgen
mapped[ri]++;
olddef = 0;
for (int j = 1; j <= pmap.Size(); j++)
olddef += sqr (nelonnode[pmap.Get(j)]);
for (auto j : pmap.Range())
olddef += sqr (nelonnode[pmap[j]]);
olddef += rule.bonus;
newdef = 0;
for (int j = 1; j <= pmap.Size(); j++)
if (rule.reused.Get(j))
newdef += sqr (nelonnode[pmap.Get(j)] +
rule.incelsonnode.Get(j));
for (auto j : pmap.Range())
if (rule.reused[j])
newdef += sqr (nelonnode[pmap[j]] +
rule.incelsonnode[j]);
if (newdef > olddef)
continue;
@ -398,7 +410,8 @@ namespace netgen
// calc metric badness
double bad1 = 0, bad2 = 0;
// SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1));
auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), &pgi.Elem(1));
// auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), &pgi.Elem(1));
auto n = geo.GetNormal(surfnr, mesh.Point(pmap.First()), pgi.Data());
for (int j = 0; j < rule.oldels.Size(); j++)
bad1 += mesh[elmap[j]].CalcJacobianBadness (mesh.Points(), n);
@ -409,7 +422,7 @@ namespace netgen
const Element2d & rnel = rule.newels.Get(j);
Element2d nel(rnel.GetNP());
for (int k = 1; k <= rnel.GetNP(); k++)
nel.PNum(k) = pmap.Get(rnel.PNum(k));
nel.PNum(k) = pmap[rnel.PNum(k)];
bad2 += nel.CalcJacobianBadness (mesh.Points(), n);
}
@ -427,8 +440,8 @@ namespace netgen
nel.SetIndex (faceindex);
for (int k = 1; k <= rnel.GetNP(); k++)
{
nel.PNum(k) = pmap.Get(rnel.PNum(k));
nel.GeomInfoPi(k) = pgi.Get(rnel.PNum(k));
nel.PNum(k) = pmap[rnel.PNum(k)];
nel.GeomInfoPi(k) = pgi[rnel.PNum(k)];
}
mesh.AddSurfaceElement(nel);
@ -437,8 +450,8 @@ namespace netgen
for (int j = 0; j < rule.oldels.Size(); j++)
mesh.Delete (elmap[j]);
for (int j = 1; j <= pmap.Size(); j++)
nelonnode[pmap.Get(j)] += rule.incelsonnode.Get(j);
for (PointIndex j : pmap.Range())
nelonnode[pmap[j]] += rule.incelsonnode[j];
used[ri]++;
}

2
libsrc/meshing/meshfunc.cpp
View File

@ -73,7 +73,7 @@ namespace netgen
m.SetLocalH(mesh.GetLocalH());
ipmap[i].SetSize(num_points);
ipmap[i] = PointIndex::INVALID;
ipmap[i] = 0; // PointIndex::INVALID;
m.SetDimension( mesh.GetDimension() );
m.SetGeometry( mesh.GetGeometry() );

18
libsrc/meshing/meshtype.hpp
View File

@ -219,10 +219,11 @@ namespace netgen
// throw Exception("illegal PointIndex, use PointIndex::INVALID instead");
#endif
}
/*
friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> ();
friend istream & operator>> (istream &, PointIndex &);
friend ostream & operator<< (ostream &, const PointIndex &);
// friend istream & operator>> (istream &, PointIndex &);
// friend ostream & operator<< (ostream &, const PointIndex &);
template <int N> friend class PointIndices;
friend PointIndex operator+ (PointIndex, int);
friend PointIndex operator+ (PointIndex, size_t);
@ -259,7 +260,6 @@ namespace netgen
void DoArchive (Archive & ar) { ar & i; }
};
/*
inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); }
inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); }
@ -284,15 +284,21 @@ namespace ngcore
namespace netgen
{
// input-output is 1-based
inline istream & operator>> (istream & ist, PointIndex & pi)
{
int i; ist >> i; pi = PointIndex(i); return ist;
// int i; ist >> i; pi = PointIndex(i); return ist;
int i; ist >> i;
// pi = PointIndex(i);
pi = IndexBASE<PointIndex>()+i-1;
return ist;
}
inline ostream & operator<< (ostream & ost, const PointIndex & pi)
{
return (ost << int(pi));
int intpi = pi - IndexBASE<PointIndex>() + 1;
return (ost << intpi);
// return (ost << int(pi));
}
template <int N> class PointIndices;