mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-24 21:10:33 +05:00
Revert "Revert "Merge branch 'parallel_meshing' into 'master'""
This reverts commit 65c5e2d244.
This commit is contained in:
Matthias Hochsteger
parent
487c32bbb4
commit
48198d232d
@ -235,13 +235,14 @@ namespace ngcore
|
||||
const function<void(TaskInfo&)> * func;
|
||||
int mynr;
|
||||
int total;
|
||||
int producing_thread;
|
||||
atomic<int> * endcnt;
|
||||
|
||||
TNestedTask () { ; }
|
||||
TNestedTask (const function<void(TaskInfo&)> & _func,
|
||||
int _mynr, int _total,
|
||||
atomic<int> & _endcnt)
|
||||
: func(&_func), mynr(_mynr), total(_total), endcnt(&_endcnt)
|
||||
atomic<int> & _endcnt, int prod_tid)
|
||||
: func(&_func), mynr(_mynr), total(_total), endcnt(&_endcnt), producing_thread(prod_tid)
|
||||
{
|
||||
;
|
||||
}
|
||||
@ -260,12 +261,14 @@ namespace ngcore
|
||||
TPToken ptoken(taskqueue);
|
||||
|
||||
int num = endcnt;
|
||||
auto tid = TaskManager::GetThreadId();
|
||||
for (int i = 0; i < num; i++)
|
||||
taskqueue.enqueue (ptoken, { afunc, i, num, endcnt });
|
||||
taskqueue.enqueue (ptoken, { afunc, i, num, endcnt, tid });
|
||||
}
|
||||
|
||||
bool TaskManager :: ProcessTask()
|
||||
{
|
||||
static Timer t("process task");
|
||||
TNestedTask task;
|
||||
TCToken ctoken(taskqueue);
|
||||
|
||||
@ -282,8 +285,14 @@ namespace ngcore
|
||||
cout << "process nested, nr = " << ti.task_nr << "/" << ti.ntasks << endl;
|
||||
}
|
||||
*/
|
||||
if(trace && task.producing_thread != ti.thread_nr)
|
||||
trace->StartTask (ti.thread_nr, t, PajeTrace::Task::ID_TIMER, task.producing_thread);
|
||||
|
||||
(*task.func)(ti);
|
||||
--*task.endcnt;
|
||||
|
||||
if(trace && task.producing_thread != ti.thread_nr)
|
||||
trace->StopTask (ti.thread_nr, t);
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
|
||||
51
libsrc/meshing/debugging.hpp
Normal file
51
libsrc/meshing/debugging.hpp
Normal file
@ -0,0 +1,51 @@
|
||||
#include <meshclass.hpp>
|
||||
|
||||
|
||||
namespace netgen
|
||||
{
|
||||
inline unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0 )
|
||||
{
|
||||
static Timer t("GetOpenElements"); RegionTimer rt(t);
|
||||
auto mesh = make_unique<Mesh>();
|
||||
*mesh = m;
|
||||
|
||||
Array<bool, PointIndex> interesting_points(mesh->GetNP());
|
||||
interesting_points = false;
|
||||
|
||||
mesh->FindOpenElements(dom);
|
||||
NgArray<Element2d> openelements;
|
||||
openelements = mesh->OpenElements();
|
||||
|
||||
for (auto & el : openelements)
|
||||
for (auto i : el.PNums())
|
||||
interesting_points[i] = true;
|
||||
|
||||
for (auto & el : mesh->VolumeElements())
|
||||
{
|
||||
int num_interesting_points = 0;
|
||||
|
||||
for (auto pi : el.PNums())
|
||||
if(interesting_points[pi])
|
||||
num_interesting_points++;
|
||||
|
||||
if(num_interesting_points==0)
|
||||
el.Delete();
|
||||
el.SetIndex(num_interesting_points);
|
||||
}
|
||||
|
||||
mesh->SetMaterial(1, "1_point");
|
||||
mesh->SetMaterial(2, "2_points");
|
||||
mesh->SetMaterial(3, "3_points");
|
||||
mesh->SetMaterial(4, "4_points");
|
||||
mesh->Compress();
|
||||
|
||||
mesh->ClearSurfaceElements();
|
||||
|
||||
for (auto & el : openelements)
|
||||
mesh->AddSurfaceElement( el );
|
||||
|
||||
return mesh;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
@ -235,9 +235,11 @@ namespace netgen
|
||||
NgArray<int> & freelist, SphereList & list,
|
||||
IndexSet & insphere, IndexSet & closesphere, Array<DelaunayTet> & newels)
|
||||
{
|
||||
static Timer t("Meshing3::AddDelaunayPoint");// RegionTimer reg(t);
|
||||
// static Timer tsearch("addpoint, search");
|
||||
// static Timer tinsert("addpoint, insert");
|
||||
static Timer t("Meshing3::AddDelaunayPoint", NoTracing, NoTiming); RegionTimer reg(t);
|
||||
static Timer tsearch("addpoint, search", NoTracing, NoTiming);
|
||||
static Timer tfind("addpoint, find all tets", NoTracing, NoTiming);
|
||||
static Timer tnewtets("addpoint, build new tets", NoTracing, NoTiming);
|
||||
static Timer tinsert("addpoint, insert", NoTracing, NoTiming);
|
||||
|
||||
/*
|
||||
find any sphere, such that newp is contained in
|
||||
@ -277,7 +279,7 @@ namespace netgen
|
||||
}
|
||||
*/
|
||||
|
||||
// tsearch.Start();
|
||||
tsearch.Start();
|
||||
double minquot{1e20};
|
||||
tettree.GetFirstIntersecting
|
||||
(newp, newp, [&](const auto pi)
|
||||
@ -300,7 +302,7 @@ namespace netgen
|
||||
}
|
||||
return false;
|
||||
} );
|
||||
// tsearch.Stop();
|
||||
tsearch.Stop();
|
||||
|
||||
if (cfelind == -1)
|
||||
{
|
||||
@ -308,6 +310,7 @@ namespace netgen
|
||||
return;
|
||||
}
|
||||
|
||||
tfind.Start();
|
||||
/*
|
||||
insphere: point is in sphere -> delete element
|
||||
closesphere: point is close to sphere -> considered for same center
|
||||
@ -399,6 +402,8 @@ namespace netgen
|
||||
}
|
||||
} // while (changed)
|
||||
|
||||
tfind.Stop();
|
||||
tnewtets.Start();
|
||||
newels.SetSize(0);
|
||||
|
||||
Element2d face(TRIG);
|
||||
@ -553,10 +558,11 @@ namespace netgen
|
||||
tpmax.SetToMax (*pp[k]);
|
||||
}
|
||||
tpmax = tpmax + 0.01 * (tpmax - tpmin);
|
||||
// tinsert.Start();
|
||||
tinsert.Start();
|
||||
tettree.Insert (tpmin, tpmax, nelind);
|
||||
// tinsert.Stop();
|
||||
tinsert.Stop();
|
||||
}
|
||||
tnewtets.Stop();
|
||||
}
|
||||
|
||||
|
||||
@ -1627,20 +1633,20 @@ namespace netgen
|
||||
// tempmesh.Save ("tempmesh.vol");
|
||||
|
||||
{
|
||||
RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0);
|
||||
for (int i = 1; i <= 4; i++)
|
||||
{
|
||||
MeshOptimize3d meshopt(mp);
|
||||
tempmesh.Compress();
|
||||
tempmesh.FindOpenElements ();
|
||||
|
||||
RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0);
|
||||
for (auto i : Range(10))
|
||||
{
|
||||
PrintMessage (5, "Num open: ", tempmesh.GetNOpenElements());
|
||||
tempmesh.CalcSurfacesOfNode ();
|
||||
|
||||
if(i%5==0)
|
||||
tempmesh.FreeOpenElementsEnvironment (1);
|
||||
|
||||
MeshOptimize3d meshopt(mp);
|
||||
// tempmesh.CalcSurfacesOfNode();
|
||||
meshopt.SwapImprove(tempmesh, OPT_CONFORM);
|
||||
}
|
||||
tempmesh.Compress();
|
||||
}
|
||||
|
||||
MeshQuality3d (tempmesh);
|
||||
|
||||
@ -1,15 +1,42 @@
|
||||
#ifndef FILE_IMPROVE2
|
||||
#define FILE_IMPROVE2
|
||||
|
||||
inline void AppendEdges( const Element2d & elem, PointIndex pi, Array<std::tuple<PointIndex,PointIndex>> & edges )
|
||||
{
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
PointIndex pi0 = elem[j];
|
||||
PointIndex pi1 = elem[(j+1)%3];
|
||||
if (pi1 < pi0) Swap(pi0, pi1);
|
||||
if(pi0==pi)
|
||||
edges.Append(std::make_tuple(pi0, pi1));
|
||||
}
|
||||
}
|
||||
|
||||
inline void AppendEdges( const Element & elem, PointIndex pi, Array<std::tuple<PointIndex,PointIndex>> & edges )
|
||||
{
|
||||
static constexpr int tetedges[6][2] =
|
||||
{ { 0, 1 }, { 0, 2 }, { 0, 3 },
|
||||
{ 1, 2 }, { 1, 3 }, { 2, 3 } };
|
||||
|
||||
if(elem.flags.fixed)
|
||||
return;
|
||||
for (int j = 0; j < 6; j++)
|
||||
{
|
||||
PointIndex pi0 = elem[tetedges[j][0]];
|
||||
PointIndex pi1 = elem[tetedges[j][1]];
|
||||
if (pi1 < pi0) Swap(pi0, pi1);
|
||||
if(pi0==pi)
|
||||
edges.Append(std::make_tuple(pi0, pi1));
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TINDEX>
|
||||
void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges )
|
||||
{
|
||||
static_assert(is_same_v<TINDEX, ElementIndex>||is_same_v<TINDEX,SurfaceElementIndex>, "Invalid type for TINDEX");
|
||||
static Timer tbuild_edges("Build edges"); RegionTimer reg(tbuild_edges);
|
||||
|
||||
static constexpr int tetedges[6][2] =
|
||||
{ { 0, 1 }, { 0, 2 }, { 0, 3 },
|
||||
{ 1, 2 }, { 1, 3 }, { 2, 3 } };
|
||||
|
||||
int ntasks = 4*ngcore::TaskManager::GetMaxThreads();
|
||||
Array<Array<std::tuple<PointIndex,PointIndex>>> task_edges(ntasks);
|
||||
|
||||
@ -26,29 +53,7 @@ void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & element
|
||||
const auto & elem = mesh[ei];
|
||||
if (elem.IsDeleted()) continue;
|
||||
|
||||
static_assert(is_same_v<TINDEX, ElementIndex>||is_same_v<TINDEX,SurfaceElementIndex>, "Invalid type for TINDEX");
|
||||
if constexpr(is_same_v<TINDEX, SurfaceElementIndex>)
|
||||
{
|
||||
for (int j = 0; j < 3; j++)
|
||||
{
|
||||
PointIndex pi0 = elem[j];
|
||||
PointIndex pi1 = elem[(j+1)%3];
|
||||
if (pi1 < pi0) Swap(pi0, pi1);
|
||||
if(pi0==pi)
|
||||
local_edges.Append(std::make_tuple(pi0, pi1));
|
||||
}
|
||||
}
|
||||
else if constexpr(is_same_v<TINDEX, ElementIndex>)
|
||||
{
|
||||
for (int j = 0; j < 6; j++)
|
||||
{
|
||||
PointIndex pi0 = elem[tetedges[j][0]];
|
||||
PointIndex pi1 = elem[tetedges[j][1]];
|
||||
if (pi1 < pi0) Swap(pi0, pi1);
|
||||
if(pi0==pi)
|
||||
local_edges.Append(std::make_tuple(pi0, pi1));
|
||||
}
|
||||
}
|
||||
AppendEdges(elem, pi, local_edges);
|
||||
}
|
||||
QuickSort(local_edges);
|
||||
|
||||
|
||||
@ -2717,8 +2717,24 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
int ne = mesh.GetNE();
|
||||
|
||||
mesh.BuildBoundaryEdges(false);
|
||||
BitArray free_points(mesh.GetNP()+PointIndex::BASE);
|
||||
free_points.Clear();
|
||||
|
||||
auto elementsonnode = mesh.CreatePoint2ElementTable();
|
||||
ParallelForRange(mesh.VolumeElements().Range(), [&] (auto myrange)
|
||||
{
|
||||
for (ElementIndex eli : myrange)
|
||||
{
|
||||
const auto & el = mesh[eli];
|
||||
if(el.flags.fixed)
|
||||
continue;
|
||||
|
||||
for (auto pi : el.PNums())
|
||||
if(!free_points[pi])
|
||||
free_points.SetBitAtomic(pi);
|
||||
}
|
||||
});
|
||||
|
||||
auto elementsonnode = mesh.CreatePoint2ElementTable(free_points);
|
||||
|
||||
NgArray<ElementIndex> hasbothpoints;
|
||||
|
||||
@ -2736,7 +2752,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
const Element2d & hel = mesh.OpenElement(i);
|
||||
INDEX_3 face(hel[0], hel[1], hel[2]);
|
||||
face.Sort();
|
||||
faces.Set (face, 1);
|
||||
faces.Set (face, i);
|
||||
}
|
||||
}
|
||||
|
||||
@ -2755,6 +2771,8 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
|
||||
tloop.Start();
|
||||
|
||||
auto num_elements_before = mesh.VolumeElements().Range().Next();
|
||||
|
||||
ParallelForRange(Range(edges), [&] (auto myrange)
|
||||
{
|
||||
for(auto i : myrange)
|
||||
@ -2786,6 +2804,30 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
|
||||
PrintMessage (5, cnt, " swaps performed");
|
||||
|
||||
if(goal == OPT_CONFORM)
|
||||
{
|
||||
// Remove open elements that were closed by new tets
|
||||
auto & open_els = mesh.OpenElements();
|
||||
|
||||
for (auto & el : mesh.VolumeElements().Range( num_elements_before, mesh.VolumeElements().Range().Next() ))
|
||||
{
|
||||
for (auto i : Range(1,5))
|
||||
{
|
||||
Element2d sel;
|
||||
el.GetFace(i, sel);
|
||||
INDEX_3 face(sel[0], sel[1], sel[2]);
|
||||
face.Sort();
|
||||
if(faces.Used(face))
|
||||
open_els[faces.Get(face)-1].Delete();
|
||||
}
|
||||
}
|
||||
|
||||
for(int i=open_els.Size()-1; i>=0; i--)
|
||||
if(open_els[i].IsDeleted())
|
||||
open_els.Delete(i);
|
||||
|
||||
mesh.DeleteBoundaryEdges();
|
||||
}
|
||||
mesh.Compress ();
|
||||
|
||||
multithread.task = savetask;
|
||||
|
||||
@ -91,6 +91,45 @@ namespace netgen
|
||||
delete root;
|
||||
}
|
||||
|
||||
unique_ptr<LocalH> LocalH :: Copy ()
|
||||
{
|
||||
static Timer t("LocalH::Copy"); RegionTimer rt(t);
|
||||
auto lh = make_unique<LocalH>(boundingbox, grading, dimension);
|
||||
std::map<GradingBox*, GradingBox*> mapping;
|
||||
lh->boxes.SetSize(boxes.Size());
|
||||
|
||||
for(auto i : boxes.Range())
|
||||
{
|
||||
lh->boxes[i] = new GradingBox();
|
||||
auto & bnew = *lh->boxes[i];
|
||||
auto & b = *boxes[i];
|
||||
bnew.xmid[0] = b.xmid[0];
|
||||
bnew.xmid[1] = b.xmid[1];
|
||||
bnew.xmid[2] = b.xmid[2];
|
||||
bnew.h2 = b.h2;
|
||||
bnew.hopt = b.hopt;
|
||||
bnew.flags = b.flags;
|
||||
mapping[&b] = &bnew;
|
||||
}
|
||||
|
||||
for(auto i : boxes.Range())
|
||||
{
|
||||
auto & bnew = *lh->boxes[i];
|
||||
auto & b = *boxes[i];
|
||||
for(auto k : Range(8))
|
||||
{
|
||||
if(b.childs[k])
|
||||
bnew.childs[k] = mapping[b.childs[k]];
|
||||
}
|
||||
|
||||
if(b.father)
|
||||
bnew.father = mapping[b.father];
|
||||
}
|
||||
|
||||
lh->root = mapping[root];
|
||||
return lh;
|
||||
}
|
||||
|
||||
void LocalH :: Delete ()
|
||||
{
|
||||
root->DeleteChilds();
|
||||
@ -405,8 +444,8 @@ namespace netgen
|
||||
void LocalH :: FindInnerBoxes (AdFront3 * adfront,
|
||||
int (*testinner)(const Point3d & p1))
|
||||
{
|
||||
static int timer = NgProfiler::CreateTimer ("LocalH::FindInnerBoxes");
|
||||
NgProfiler::RegionTimer reg (timer);
|
||||
static Timer timer("LocalH::FindInnerBoxes");
|
||||
RegionTimer reg (timer);
|
||||
|
||||
|
||||
int nf = adfront->GetNF();
|
||||
@ -812,6 +851,7 @@ namespace netgen
|
||||
|
||||
void LocalH :: GetOuterPoints (NgArray<Point<3> > & points)
|
||||
{
|
||||
static Timer t("LocalH::GetOuterPoints"); RegionTimer rt(t);
|
||||
for (int i = 0; i < boxes.Size(); i++)
|
||||
if (!boxes[i]->flags.isinner && !boxes[i]->flags.cutboundary)
|
||||
points.Append ( boxes[i] -> PMid());
|
||||
|
||||
@ -98,6 +98,8 @@ namespace netgen
|
||||
|
||||
~LocalH();
|
||||
///
|
||||
unique_ptr<LocalH> Copy();
|
||||
///
|
||||
void Delete();
|
||||
///
|
||||
void DoArchive(Archive& ar);
|
||||
|
||||
@ -384,6 +384,8 @@ namespace netgen
|
||||
volelements.Append (el);
|
||||
}
|
||||
volelements.Last().flags.illegal_valid = 0;
|
||||
volelements.Last().flags.fixed = 0;
|
||||
volelements.Last().flags.deleted = 0;
|
||||
|
||||
// while (volelements.Size() > eltyps.Size())
|
||||
// eltyps.Append (FREEELEMENT);
|
||||
@ -405,6 +407,8 @@ namespace netgen
|
||||
|
||||
volelements[ei] = el;
|
||||
volelements[ei].flags.illegal_valid = 0;
|
||||
volelements[ei].flags.fixed = 0;
|
||||
volelements[ei].flags.deleted = 0;
|
||||
}
|
||||
|
||||
|
||||
@ -436,6 +440,7 @@ namespace netgen
|
||||
|
||||
void Mesh :: Save (ostream & outfile) const
|
||||
{
|
||||
static Timer timer("Mesh::Save"); RegionTimer rt(timer);
|
||||
int i, j;
|
||||
|
||||
double scale = 1; // globflags.GetNumFlag ("scale", 1);
|
||||
@ -2910,6 +2915,7 @@ namespace netgen
|
||||
|
||||
void Mesh :: FreeOpenElementsEnvironment (int layers)
|
||||
{
|
||||
static Timer timer("FreeOpenElementsEnvironment"); RegionTimer rt(timer);
|
||||
int i, j, k;
|
||||
PointIndex pi;
|
||||
const int large = 9999;
|
||||
@ -6570,6 +6576,8 @@ namespace netgen
|
||||
[&](auto & table, ElementIndex ei)
|
||||
{
|
||||
const auto & el = (*this)[ei];
|
||||
if(el.IsDeleted())
|
||||
return;
|
||||
|
||||
for (PointIndex pi : el.PNums())
|
||||
if(free_points[pi])
|
||||
@ -6581,6 +6589,8 @@ namespace netgen
|
||||
[&](auto & table, ElementIndex ei)
|
||||
{
|
||||
const auto & el = (*this)[ei];
|
||||
if(el.IsDeleted())
|
||||
return;
|
||||
|
||||
for (PointIndex pi : el.PNums())
|
||||
table.Add (pi, ei);
|
||||
|
||||
@ -487,6 +487,8 @@ namespace netgen
|
||||
|
||||
auto & OpenElements() const { return openelements; }
|
||||
|
||||
auto & OpenElements() { return openelements; }
|
||||
|
||||
/// are also quads open elements
|
||||
bool HasOpenQuads () const;
|
||||
|
||||
@ -510,6 +512,11 @@ namespace netgen
|
||||
return boundaryedges->Used (i2);
|
||||
}
|
||||
|
||||
void DeleteBoundaryEdges ()
|
||||
{
|
||||
boundaryedges = nullptr;
|
||||
}
|
||||
|
||||
bool IsSegment (PointIndex pi1, PointIndex pi2) const
|
||||
{
|
||||
INDEX_2 i2 (pi1, pi2);
|
||||
|
||||
@ -1,5 +1,6 @@
|
||||
#include <mystdlib.h>
|
||||
#include "meshing.hpp"
|
||||
#include "debugging.hpp"
|
||||
|
||||
namespace netgen
|
||||
{
|
||||
@ -10,80 +11,405 @@ namespace netgen
|
||||
extern const char * pyramidrules2[];
|
||||
extern const char * hexrules[];
|
||||
|
||||
|
||||
// extern double teterrpow;
|
||||
MESHING3_RESULT MeshVolume (const MeshingParameters & c_mp, Mesh& mesh3d)
|
||||
struct MeshingData
|
||||
{
|
||||
static Timer t("MeshVolume"); RegionTimer reg(t);
|
||||
int domain;
|
||||
|
||||
MeshingParameters mp = c_mp; // copy mp to change them here
|
||||
int oldne;
|
||||
int meshed;
|
||||
// mesh for one domain (contains all adjacent surface elments)
|
||||
unique_ptr<Mesh> mesh;
|
||||
|
||||
// maps from local (domain) mesh to global mesh
|
||||
Array<PointIndex, PointIndex> pmap;
|
||||
|
||||
Array<INDEX_2> connected_pairs;
|
||||
|
||||
MeshingParameters mp;
|
||||
|
||||
unique_ptr<Meshing3> meshing;
|
||||
};
|
||||
|
||||
// extract surface meshes belonging to individual domains
|
||||
Array<MeshingData> DivideMesh(Mesh & mesh, const MeshingParameters & mp)
|
||||
{
|
||||
static Timer timer("DivideMesh"); RegionTimer rt(timer);
|
||||
|
||||
Array<MeshingData> ret;
|
||||
auto num_domains = mesh.GetNDomains();
|
||||
|
||||
ret.SetSize(num_domains);
|
||||
if(num_domains==1 || mp.only3D_domain_nr)
|
||||
{
|
||||
// no need to divide mesh, just fill in meshing data
|
||||
ret[0].domain = 1;
|
||||
if(mp.only3D_domain_nr)
|
||||
ret[0].domain = mp.only3D_domain_nr;
|
||||
|
||||
ret[0].mesh.reset(&mesh); // careful, this unique_ptr must not delete &mesh! (it will be released in MergeMeshes after meshing)
|
||||
ret[0].mp = mp;
|
||||
return ret;
|
||||
}
|
||||
|
||||
Array<Array<PointIndex, PointIndex>> ipmap;
|
||||
ipmap.SetSize(num_domains);
|
||||
auto dim = mesh.GetDimension();
|
||||
auto num_points = mesh.GetNP();
|
||||
auto num_facedescriptors = mesh.GetNFD();
|
||||
|
||||
auto & identifications = mesh.GetIdentifications();
|
||||
|
||||
for(auto i : Range(ret))
|
||||
{
|
||||
auto & md = ret[i];
|
||||
md.domain = i+1;
|
||||
|
||||
md.mp = mp;
|
||||
md.mp.maxh = min2 (mp.maxh, mesh.MaxHDomain(md.domain));
|
||||
|
||||
ret[i].mesh = make_unique<Mesh>();
|
||||
auto & m = *ret[i].mesh;
|
||||
|
||||
m.SetLocalH(mesh.GetLocalH());
|
||||
|
||||
ipmap[i].SetSize(num_points);
|
||||
ipmap[i] = PointIndex::INVALID;
|
||||
m.SetDimension( mesh.GetDimension() );
|
||||
|
||||
for(auto i : Range(1, num_facedescriptors+1))
|
||||
m.AddFaceDescriptor( mesh.GetFaceDescriptor(i) );
|
||||
}
|
||||
|
||||
// mark used points for each domain, add surface elements (with wrong point numbers) to domain mesh
|
||||
for(const auto & sel : mesh.SurfaceElements())
|
||||
{
|
||||
const auto & fd = mesh.GetFaceDescriptor(sel.GetIndex());
|
||||
int dom_in = fd.DomainIn();
|
||||
int dom_out = fd.DomainOut();
|
||||
|
||||
for( auto dom : {dom_in, dom_out} )
|
||||
{
|
||||
if(dom==0)
|
||||
continue;
|
||||
|
||||
auto & sels = ret[dom-1].mesh->SurfaceElements();
|
||||
for(auto pi : sel.PNums())
|
||||
ipmap[dom-1][pi] = 1;
|
||||
sels.Append(sel);
|
||||
}
|
||||
}
|
||||
|
||||
// mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points?
|
||||
for(auto pi : mesh.LockedPoints())
|
||||
for(auto i : Range(ret))
|
||||
ipmap[i][pi] = 1;
|
||||
|
||||
// add used points to domain mesh, build point mapping
|
||||
for(auto i : Range(ret))
|
||||
{
|
||||
auto & m = *ret[i].mesh;
|
||||
auto & pmap = ret[i].pmap;
|
||||
|
||||
for(auto pi : Range(ipmap[i]))
|
||||
if(ipmap[i][pi])
|
||||
{
|
||||
auto pi_new = m.AddPoint( mesh[pi] );
|
||||
ipmap[i][pi] = pi_new;
|
||||
pmap.Append( pi );
|
||||
}
|
||||
}
|
||||
|
||||
NgArray<INDEX_2> connectednodes;
|
||||
|
||||
if (!mesh3d.HasLocalHFunction()) mesh3d.CalcLocalH(mp.grading);
|
||||
|
||||
mesh3d.Compress();
|
||||
|
||||
// mesh3d.PrintMemInfo (cout);
|
||||
|
||||
if (mp.checkoverlappingboundary)
|
||||
if (mesh3d.CheckOverlappingBoundary())
|
||||
throw NgException ("Stop meshing since boundary mesh is overlapping");
|
||||
|
||||
int nonconsist = 0;
|
||||
for (int k = 1; k <= mesh3d.GetNDomains(); k++)
|
||||
for(auto i : Range(ret))
|
||||
{
|
||||
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
|
||||
continue;
|
||||
PrintMessage (3, "Check subdomain ", k, " / ", mesh3d.GetNDomains());
|
||||
auto & imap = ipmap[i];
|
||||
auto & m = *ret[i].mesh;
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
for (auto & sel : m.SurfaceElements())
|
||||
for(auto & pi : sel.PNums())
|
||||
pi = imap[pi];
|
||||
|
||||
/*
|
||||
bool res = mesh3d.CheckOverlappingBoundary();
|
||||
if (res)
|
||||
for (int nr = 1; nr <= identifications.GetMaxNr(); nr++)
|
||||
if (identifications.GetType(nr) != Identifications::PERIODIC)
|
||||
{
|
||||
PrintError ("Surface is overlapping !!");
|
||||
nonconsist = 1;
|
||||
identifications.GetPairs (nr, connectednodes);
|
||||
for (auto pair : connectednodes)
|
||||
{
|
||||
auto pi0 = pair[0];
|
||||
auto pi1 = pair[1];
|
||||
if(imap[pi0].IsValid() && imap[pi1].IsValid())
|
||||
ret[i].connected_pairs.Append({imap[pi0], imap[pi1]});
|
||||
}
|
||||
}
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
*/
|
||||
|
||||
bool res = (mesh3d.CheckConsistentBoundary() != 0);
|
||||
if (res)
|
||||
void CloseOpenQuads( MeshingData & md)
|
||||
{
|
||||
PrintError ("Surface mesh not consistent");
|
||||
nonconsist = 1;
|
||||
auto & mesh = *md.mesh;
|
||||
auto domain = md.domain;
|
||||
MeshingParameters & mp = md.mp;
|
||||
|
||||
int oldne;
|
||||
if (multithread.terminate)
|
||||
return;
|
||||
|
||||
mesh.CalcSurfacesOfNode();
|
||||
mesh.FindOpenElements(domain);
|
||||
|
||||
if (!mesh.GetNOpenElements())
|
||||
return;
|
||||
|
||||
for (int qstep = 0; qstep <= 3; qstep++)
|
||||
{
|
||||
if (qstep == 0 && !mp.try_hexes) continue;
|
||||
|
||||
if (mesh.HasOpenQuads())
|
||||
{
|
||||
string rulefile = ngdir;
|
||||
|
||||
const char ** rulep = NULL;
|
||||
switch (qstep)
|
||||
{
|
||||
case 0:
|
||||
rulep = hexrules;
|
||||
break;
|
||||
case 1:
|
||||
rulep = prismrules2;
|
||||
break;
|
||||
case 2: // connect pyramid to triangle
|
||||
rulep = pyramidrules2;
|
||||
break;
|
||||
case 3: // connect to vis-a-vis point
|
||||
rulep = pyramidrules;
|
||||
break;
|
||||
}
|
||||
|
||||
Meshing3 meshing(rulep);
|
||||
|
||||
MeshingParameters mpquad = mp;
|
||||
|
||||
mpquad.giveuptol = 15;
|
||||
mpquad.baseelnp = 4;
|
||||
mpquad.starshapeclass = 1000;
|
||||
mpquad.check_impossible = qstep == 1; // for prisms only (air domain in trafo)
|
||||
|
||||
|
||||
for (PointIndex pi : mesh.Points().Range())
|
||||
meshing.AddPoint (mesh[pi], pi);
|
||||
|
||||
for (auto pair : md.connected_pairs)
|
||||
meshing.AddConnectedPair (pair);
|
||||
|
||||
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
|
||||
{
|
||||
Element2d hel = mesh.OpenElement(i);
|
||||
meshing.AddBoundaryElement (hel);
|
||||
}
|
||||
|
||||
oldne = mesh.GetNE();
|
||||
|
||||
meshing.GenerateMesh (mesh, mpquad);
|
||||
|
||||
for (int i = oldne + 1; i <= mesh.GetNE(); i++)
|
||||
mesh.VolumeElement(i).SetIndex (domain);
|
||||
|
||||
(*testout)
|
||||
<< "mesh has " << mesh.GetNE() << " prism/pyramid elements" << endl;
|
||||
|
||||
mesh.FindOpenElements(domain);
|
||||
}
|
||||
}
|
||||
|
||||
if (nonconsist)
|
||||
|
||||
if (mesh.HasOpenQuads())
|
||||
{
|
||||
PrintError ("Stop meshing since surface mesh not consistent");
|
||||
throw NgException ("Stop meshing since surface mesh not consistent");
|
||||
PrintSysError ("mesh has still open quads");
|
||||
throw NgException ("Stop meshing since too many attempts");
|
||||
// return MESHING3_GIVEUP;
|
||||
}
|
||||
}
|
||||
|
||||
double globmaxh = mp.maxh;
|
||||
|
||||
for (int k = 1; k <= mesh3d.GetNDomains(); k++)
|
||||
void PrepareForBlockFillLocalH(MeshingData & md)
|
||||
{
|
||||
static Timer t("PrepareForBlockFillLocalH"); RegionTimer rt(t);
|
||||
md.meshing = make_unique<Meshing3>(nullptr);
|
||||
|
||||
auto & mesh = *md.mesh;
|
||||
|
||||
mesh.CalcSurfacesOfNode();
|
||||
mesh.FindOpenElements(md.domain);
|
||||
|
||||
for (PointIndex pi : mesh.Points().Range())
|
||||
md.meshing->AddPoint (mesh[pi], pi);
|
||||
|
||||
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
|
||||
md.meshing->AddBoundaryElement (mesh.OpenElement(i));
|
||||
|
||||
if (mesh.HasLocalHFunction())
|
||||
md.meshing->PrepareBlockFillLocalH(mesh, md.mp);
|
||||
}
|
||||
|
||||
|
||||
void MeshDomain( MeshingData & md)
|
||||
{
|
||||
auto & mesh = *md.mesh;
|
||||
auto domain = md.domain;
|
||||
MeshingParameters & mp = md.mp;
|
||||
|
||||
if (mp.delaunay && mesh.GetNOpenElements())
|
||||
{
|
||||
int oldne = mesh.GetNE();
|
||||
|
||||
md.meshing->Delaunay (mesh, domain, mp);
|
||||
|
||||
for (int i = oldne + 1; i <= mesh.GetNE(); i++)
|
||||
mesh.VolumeElement(i).SetIndex (domain);
|
||||
|
||||
PrintMessage (3, mesh.GetNP(), " points, ",
|
||||
mesh.GetNE(), " elements");
|
||||
}
|
||||
|
||||
Box<3> domain_bbox( Box<3>::EMPTY_BOX );
|
||||
|
||||
for (auto & sel : mesh.SurfaceElements())
|
||||
{
|
||||
if (sel.IsDeleted() ) continue;
|
||||
|
||||
for (auto pi : sel.PNums())
|
||||
domain_bbox.Add (mesh[pi]);
|
||||
}
|
||||
domain_bbox.Increase (0.01 * domain_bbox.Diam());
|
||||
|
||||
mesh.FindOpenElements(domain);
|
||||
|
||||
int cntsteps = 0;
|
||||
int meshed;
|
||||
if (mesh.GetNOpenElements())
|
||||
do
|
||||
{
|
||||
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
|
||||
continue;
|
||||
if (multithread.terminate)
|
||||
break;
|
||||
|
||||
mesh.FindOpenElements(domain);
|
||||
PrintMessage (5, mesh.GetNOpenElements(), " open faces");
|
||||
GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol");
|
||||
cntsteps++;
|
||||
|
||||
|
||||
if (cntsteps > mp.maxoutersteps)
|
||||
throw NgException ("Stop meshing since too many attempts");
|
||||
|
||||
PrintMessage (1, "start tetmeshing");
|
||||
|
||||
Meshing3 meshing(tetrules);
|
||||
|
||||
Array<PointIndex, PointIndex> glob2loc(mesh.GetNP());
|
||||
glob2loc = PointIndex::INVALID;
|
||||
|
||||
for (PointIndex pi : mesh.Points().Range())
|
||||
if (domain_bbox.IsIn (mesh[pi]))
|
||||
glob2loc[pi] = meshing.AddPoint (mesh[pi], pi);
|
||||
|
||||
for (auto sel : mesh.OpenElements() )
|
||||
{
|
||||
for(auto & pi : sel.PNums())
|
||||
pi = glob2loc[pi];
|
||||
meshing.AddBoundaryElement (sel);
|
||||
}
|
||||
|
||||
int oldne = mesh.GetNE();
|
||||
|
||||
mp.giveuptol = 15 + 10 * cntsteps;
|
||||
mp.sloppy = 5;
|
||||
meshing.GenerateMesh (mesh, mp);
|
||||
|
||||
for (ElementIndex ei = oldne; ei < mesh.GetNE(); ei++)
|
||||
mesh[ei].SetIndex (domain);
|
||||
|
||||
|
||||
mesh.CalcSurfacesOfNode();
|
||||
mesh.FindOpenElements(domain);
|
||||
|
||||
// teterrpow = 2;
|
||||
if (mesh.GetNOpenElements() != 0)
|
||||
{
|
||||
meshed = 0;
|
||||
PrintMessage (5, mesh.GetNOpenElements(), " open faces found");
|
||||
|
||||
MeshOptimize3d optmesh(mp);
|
||||
|
||||
const char * optstr = "mcmstmcmstmcmstmcm";
|
||||
for (size_t j = 1; j <= strlen(optstr); j++)
|
||||
{
|
||||
mesh.CalcSurfacesOfNode();
|
||||
mesh.FreeOpenElementsEnvironment(2);
|
||||
mesh.CalcSurfacesOfNode();
|
||||
|
||||
switch (optstr[j-1])
|
||||
{
|
||||
case 'c': optmesh.CombineImprove(mesh, OPT_REST); break;
|
||||
case 'd': optmesh.SplitImprove(mesh, OPT_REST); break;
|
||||
case 's': optmesh.SwapImprove(mesh, OPT_REST); break;
|
||||
case 't': optmesh.SwapImprove2(mesh, OPT_REST); break;
|
||||
case 'm': mesh.ImproveMesh(mp, OPT_REST); break;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
mesh.FindOpenElements();
|
||||
PrintMessage (3, "Call remove problem");
|
||||
// mesh.Save("before_remove.vol");
|
||||
RemoveProblem (mesh, domain);
|
||||
// mesh.Save("after_remove.vol");
|
||||
mesh.FindOpenElements();
|
||||
}
|
||||
else
|
||||
{
|
||||
meshed = 1;
|
||||
PrintMessage (1, "Success !");
|
||||
}
|
||||
}
|
||||
while (!meshed);
|
||||
|
||||
{
|
||||
PrintMessage (3, "Check subdomain ", domain, " / ", mesh.GetNDomains());
|
||||
|
||||
mesh.FindOpenElements(domain);
|
||||
|
||||
bool res = (mesh.CheckConsistentBoundary() != 0);
|
||||
if (res)
|
||||
{
|
||||
// mesh.Save("output.vol");
|
||||
PrintError ("Surface mesh not consistent");
|
||||
throw NgException ("Stop meshing since surface mesh not consistent");
|
||||
}
|
||||
}
|
||||
// OptimizeVolume( md.mp, mesh );
|
||||
}
|
||||
|
||||
void MeshDomain(Mesh & mesh3d, const MeshingParameters & c_mp, int k, const Identifications & identifications)
|
||||
{
|
||||
MeshingParameters mp = c_mp; // copy mp to change them here
|
||||
NgArray<INDEX_2> connectednodes;
|
||||
|
||||
int oldne;
|
||||
int meshed;
|
||||
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
|
||||
return;
|
||||
if (multithread.terminate)
|
||||
return;
|
||||
|
||||
PrintMessage (2, "");
|
||||
PrintMessage (1, "Meshing subdomain ", k, " of ", mesh3d.GetNDomains());
|
||||
(*testout) << "Meshing subdomain " << k << endl;
|
||||
|
||||
mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
|
||||
mp.maxh = min2 (mp.maxh, mesh3d.MaxHDomain(k));
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
if (!mesh3d.GetNOpenElements())
|
||||
continue;
|
||||
return;
|
||||
|
||||
|
||||
|
||||
@ -154,13 +480,13 @@ namespace netgen
|
||||
for (int i = 1; i <= connectednodes.Size(); i++)
|
||||
meshing.AddConnectedPair (connectednodes.Get(i));
|
||||
*/
|
||||
for (int nr = 1; nr <= mesh3d.GetIdentifications().GetMaxNr(); nr++)
|
||||
if (mesh3d.GetIdentifications().GetType(nr) != Identifications::PERIODIC)
|
||||
{
|
||||
mesh3d.GetIdentifications().GetPairs (nr, connectednodes);
|
||||
for (auto pair : connectednodes)
|
||||
meshing.AddConnectedPair (pair);
|
||||
}
|
||||
// for (int nr = 1; nr <= identifications.GetMaxNr(); nr++)
|
||||
// if (identifications.GetType(nr) != Identifications::PERIODIC)
|
||||
// {
|
||||
// identifications.GetPairs (nr, connectednodes);
|
||||
// for (auto pair : connectednodes)
|
||||
// meshing.AddConnectedPair (pair);
|
||||
// }
|
||||
|
||||
for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
{
|
||||
@ -311,334 +637,122 @@ namespace netgen
|
||||
|
||||
PrintMessage (1, mesh3d.GetNP(), " points, ",
|
||||
mesh3d.GetNE(), " elements");
|
||||
}
|
||||
|
||||
mp.maxh = globmaxh;
|
||||
|
||||
MeshQuality3d (mesh3d);
|
||||
|
||||
return MESHING3_OK;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/*
|
||||
|
||||
|
||||
MESHING3_RESULT MeshVolumeOld (MeshingParameters & mp, Mesh& mesh3d)
|
||||
{
|
||||
int i, k, oldne;
|
||||
|
||||
|
||||
int meshed;
|
||||
int cntsteps;
|
||||
|
||||
|
||||
PlotStatistics3d * pstat;
|
||||
if (globflags.GetNumFlag("silentflag", 1) <= 2)
|
||||
pstat = new XPlotStatistics3d;
|
||||
else
|
||||
pstat = new TerminalPlotStatistics3d;
|
||||
|
||||
cntsteps = 0;
|
||||
do
|
||||
{
|
||||
cntsteps++;
|
||||
if (cntsteps > mp.maxoutersteps)
|
||||
{
|
||||
return MESHING3_OUTERSTEPSEXCEEDED;
|
||||
}
|
||||
|
||||
|
||||
int noldp = mesh3d.GetNP();
|
||||
|
||||
|
||||
if ( (cntsteps == 1) && globflags.GetDefineFlag ("delaunay"))
|
||||
{
|
||||
cntsteps ++;
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
|
||||
|
||||
for (k = 1; k <= mesh3d.GetNDomains(); k++)
|
||||
{
|
||||
Meshing3 meshing(NULL, pstat);
|
||||
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
|
||||
return;
|
||||
PrintMessage (3, "Check subdomain ", k, " / ", mesh3d.GetNDomains());
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
for (i = 1; i <= noldp; i++)
|
||||
meshing.AddPoint (mesh3d.Point(i), i);
|
||||
|
||||
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
bool res = (mesh3d.CheckConsistentBoundary() != 0);
|
||||
if (res)
|
||||
{
|
||||
if (mesh3d.OpenElement(i).GetIndex() == k)
|
||||
meshing.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
PrintError ("Surface mesh not consistent");
|
||||
throw NgException ("Stop meshing since surface mesh not consistent");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
oldne = mesh3d.GetNE();
|
||||
if (globflags.GetDefineFlag ("blockfill"))
|
||||
void MergeMeshes( Mesh & mesh, Array<MeshingData> & md )
|
||||
{
|
||||
if (!globflags.GetDefineFlag ("localh"))
|
||||
meshing.BlockFill
|
||||
(mesh3d, mp.h * globflags.GetNumFlag ("relblockfillh", 1));
|
||||
else
|
||||
meshing.BlockFillLocalH (mesh3d);
|
||||
}
|
||||
|
||||
MeshingParameters mpd;
|
||||
meshing.Delaunay (mesh3d, mpd);
|
||||
|
||||
for (i = oldne + 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (k);
|
||||
}
|
||||
}
|
||||
|
||||
noldp = mesh3d.GetNP();
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements();
|
||||
for (k = 1; k <= mesh3d.GetNDomains(); k++)
|
||||
// todo: optimize: count elements, alloc all memory, copy vol elements in parallel
|
||||
static Timer t("MergeMeshes"); RegionTimer rt(t);
|
||||
if(md.Size()==1)
|
||||
{
|
||||
Meshing3 meshing(globflags.GetStringFlag ("rules3d", NULL), pstat);
|
||||
// assume that mesh was never divided, no need to do anything
|
||||
if(&mesh != md[0].mesh.get())
|
||||
throw Exception("Illegal Mesh pointer in MeshingData");
|
||||
|
||||
Point3d pmin, pmax;
|
||||
mesh3d.GetBox (pmin, pmax, k);
|
||||
md[0].mesh.release();
|
||||
return;
|
||||
}
|
||||
|
||||
rot.SetCenter (Center (pmin, pmax));
|
||||
|
||||
for (i = 1; i <= noldp; i++)
|
||||
meshing.AddPoint (mesh3d.Point(i), i);
|
||||
|
||||
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
for(auto & m_ : md)
|
||||
{
|
||||
if (mesh3d.OpenElement(i).GetIndex() == k)
|
||||
meshing.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
}
|
||||
auto first_new_pi = m_.pmap.Range().Next();
|
||||
auto & m = *m_.mesh;
|
||||
Array<PointIndex, PointIndex> pmap(m.Points().Size());
|
||||
for(auto pi : Range(PointIndex(PointIndex::BASE), first_new_pi))
|
||||
pmap[pi] = m_.pmap[pi];
|
||||
|
||||
oldne = mesh3d.GetNE();
|
||||
for (auto pi : Range(first_new_pi, m.Points().Range().Next()))
|
||||
pmap[pi] = mesh.AddPoint(m[pi]);
|
||||
|
||||
|
||||
if ( (cntsteps == 1) && globflags.GetDefineFlag ("blockfill"))
|
||||
for ( auto el : m.VolumeElements() )
|
||||
{
|
||||
if (!globflags.GetDefineFlag ("localh"))
|
||||
for (auto i : Range(el.GetNP()))
|
||||
el[i] = pmap[el[i]];
|
||||
el.SetIndex(m_.domain);
|
||||
mesh.AddVolumeElement(el);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void MergeMeshes( Mesh & mesh, FlatArray<Mesh> meshes, PointIndex first_new_pi )
|
||||
{
|
||||
meshing.BlockFill
|
||||
(mesh3d,
|
||||
mp.h * globflags.GetNumFlag ("relblockfillh", 1));
|
||||
}
|
||||
else
|
||||
// todo: optimize: count elements, alloc all memory, copy vol elements in parallel
|
||||
static Timer t("MergeMeshes"); RegionTimer rt(t);
|
||||
for(auto & m : meshes)
|
||||
{
|
||||
meshing.BlockFillLocalH (mesh3d);
|
||||
}
|
||||
}
|
||||
Array<PointIndex, PointIndex> pmap(m.Points().Size());
|
||||
for(auto pi : Range(PointIndex(PointIndex::BASE), first_new_pi))
|
||||
pmap[pi] = pi;
|
||||
|
||||
for (auto pi : Range(first_new_pi, m.Points().Range().Next()))
|
||||
pmap[pi] = mesh.AddPoint(m[pi]);
|
||||
|
||||
|
||||
mp.giveuptol = int(globflags.GetNumFlag ("giveuptol", 15));
|
||||
|
||||
meshing.GenerateMesh (mesh3d, mp);
|
||||
|
||||
for (i = oldne + 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (k);
|
||||
}
|
||||
|
||||
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements();
|
||||
|
||||
teterrpow = 2;
|
||||
if (mesh3d.GetNOpenElements() != 0)
|
||||
for ( auto el : m.VolumeElements() )
|
||||
{
|
||||
meshed = 0;
|
||||
(*mycout) << "Open elements found, old" << endl;
|
||||
const char * optstr = "mcmcmcmcm";
|
||||
int j;
|
||||
for (j = 1; j <= strlen(optstr); j++)
|
||||
switch (optstr[j-1])
|
||||
{
|
||||
case 'c': mesh3d.CombineImprove(); break;
|
||||
case 'd': mesh3d.SplitImprove(); break;
|
||||
case 's': mesh3d.SwapImprove(); break;
|
||||
case 'm': mesh3d.ImproveMesh(2); break;
|
||||
for (auto i : Range(el.GetNP()))
|
||||
el[i] = pmap[el[i]];
|
||||
mesh.AddVolumeElement(el);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
(*mycout) << "Call remove" << endl;
|
||||
RemoveProblem (mesh3d);
|
||||
(*mycout) << "Problem removed" << endl;
|
||||
}
|
||||
else
|
||||
meshed = 1;
|
||||
}
|
||||
while (!meshed);
|
||||
// extern double teterrpow;
|
||||
MESHING3_RESULT MeshVolume (const MeshingParameters & mp, Mesh& mesh3d)
|
||||
{
|
||||
static Timer t("MeshVolume"); RegionTimer reg(t);
|
||||
|
||||
mesh3d.Compress();
|
||||
|
||||
if (mp.checkoverlappingboundary)
|
||||
if (mesh3d.CheckOverlappingBoundary())
|
||||
throw NgException ("Stop meshing since boundary mesh is overlapping");
|
||||
|
||||
|
||||
if(mesh3d.GetNDomains()==0)
|
||||
return MESHING3_OK;
|
||||
|
||||
// localh function is built for each domain separately in blockfill ( more efficient )
|
||||
if (!mesh3d.HasLocalHFunction() && !mp.blockfill)
|
||||
mesh3d.CalcLocalH(mp.grading);
|
||||
|
||||
auto md = DivideMesh(mesh3d, mp);
|
||||
|
||||
ParallelFor( md.Range(), [&](int i)
|
||||
{
|
||||
CloseOpenQuads( md[i] );
|
||||
});
|
||||
|
||||
for(auto & md_ : md)
|
||||
PrepareForBlockFillLocalH(md_);
|
||||
|
||||
ParallelFor( md.Range(), [&](int i)
|
||||
{
|
||||
MeshDomain(md[i]);
|
||||
});
|
||||
|
||||
MergeMeshes(mesh3d, md);
|
||||
|
||||
MeshQuality3d (mesh3d);
|
||||
|
||||
return MESHING3_OK;
|
||||
}
|
||||
|
||||
*/
|
||||
|
||||
|
||||
|
||||
|
||||
/*
|
||||
MESHING3_RESULT MeshMixedVolume(MeshingParameters & mp, Mesh& mesh3d)
|
||||
{
|
||||
int i, j;
|
||||
MESHING3_RESULT res;
|
||||
Point3d pmin, pmax;
|
||||
|
||||
mp.giveuptol = 10;
|
||||
mp.baseelnp = 4;
|
||||
mp.starshapeclass = 100;
|
||||
|
||||
// TerminalPlotStatistics3d pstat;
|
||||
|
||||
Meshing3 meshing1("pyramids.rls");
|
||||
for (i = 1; i <= mesh3d.GetNP(); i++)
|
||||
meshing1.AddPoint (mesh3d.Point(i), i);
|
||||
|
||||
mesh3d.FindOpenElements();
|
||||
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
if (mesh3d.OpenElement(i).GetIndex() == 1)
|
||||
meshing1.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
|
||||
res = meshing1.GenerateMesh (mesh3d, mp);
|
||||
|
||||
mesh3d.GetBox (pmin, pmax);
|
||||
PrintMessage (1, "Mesh pyramids, res = ", res);
|
||||
if (res)
|
||||
exit (1);
|
||||
|
||||
|
||||
for (i = 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (1);
|
||||
|
||||
// do delaunay
|
||||
|
||||
mp.baseelnp = 0;
|
||||
mp.starshapeclass = 5;
|
||||
|
||||
Meshing3 meshing2(NULL);
|
||||
for (i = 1; i <= mesh3d.GetNP(); i++)
|
||||
meshing2.AddPoint (mesh3d.Point(i), i);
|
||||
|
||||
mesh3d.FindOpenElements();
|
||||
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
if (mesh3d.OpenElement(i).GetIndex() == 1)
|
||||
meshing2.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
|
||||
MeshingParameters mpd;
|
||||
meshing2.Delaunay (mesh3d, mpd);
|
||||
|
||||
for (i = 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (1);
|
||||
|
||||
|
||||
mp.baseelnp = 0;
|
||||
mp.giveuptol = 10;
|
||||
|
||||
for (int trials = 1; trials <= 50; trials++)
|
||||
{
|
||||
if (multithread.terminate)
|
||||
return MESHING3_TERMINATE;
|
||||
|
||||
Meshing3 meshing3("tetra.rls");
|
||||
for (i = 1; i <= mesh3d.GetNP(); i++)
|
||||
meshing3.AddPoint (mesh3d.Point(i), i);
|
||||
|
||||
mesh3d.FindOpenElements();
|
||||
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
if (mesh3d.OpenElement(i).GetIndex() == 1)
|
||||
meshing3.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
|
||||
if (trials > 1)
|
||||
CheckSurfaceMesh2 (mesh3d);
|
||||
res = meshing3.GenerateMesh (mesh3d, mp);
|
||||
|
||||
for (i = 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (1);
|
||||
|
||||
if (res == 0) break;
|
||||
|
||||
|
||||
|
||||
for (i = 1; i <= mesh3d.GetNE(); i++)
|
||||
{
|
||||
const Element & el = mesh3d.VolumeElement(i);
|
||||
if (el.GetNP() != 4)
|
||||
{
|
||||
for (j = 1; j <= el.GetNP(); j++)
|
||||
mesh3d.AddLockedPoint (el.PNum(j));
|
||||
}
|
||||
}
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements();
|
||||
|
||||
MeshOptimize3d optmesh;
|
||||
|
||||
teterrpow = 2;
|
||||
const char * optstr = "mcmcmcmcm";
|
||||
for (j = 1; j <= strlen(optstr); j++)
|
||||
switch (optstr[j-1])
|
||||
{
|
||||
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
|
||||
case 'd': optmesh.SplitImprove(mesh3d); break;
|
||||
case 's': optmesh.SwapImprove(mesh3d); break;
|
||||
case 'm': mesh3d.ImproveMesh(); break;
|
||||
}
|
||||
|
||||
RemoveProblem (mesh3d);
|
||||
}
|
||||
|
||||
|
||||
PrintMessage (1, "Meshing tets, res = ", res);
|
||||
if (res)
|
||||
{
|
||||
mesh3d.FindOpenElements();
|
||||
PrintSysError (1, "Open elements: ", mesh3d.GetNOpenElements());
|
||||
exit (1);
|
||||
}
|
||||
|
||||
|
||||
|
||||
for (i = 1; i <= mesh3d.GetNE(); i++)
|
||||
{
|
||||
const Element & el = mesh3d.VolumeElement(i);
|
||||
if (el.GetNP() != 4)
|
||||
{
|
||||
for (j = 1; j <= el.GetNP(); j++)
|
||||
mesh3d.AddLockedPoint (el.PNum(j));
|
||||
}
|
||||
}
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements();
|
||||
|
||||
MeshOptimize3d optmesh;
|
||||
|
||||
teterrpow = 2;
|
||||
const char * optstr = "mcmcmcmcm";
|
||||
for (j = 1; j <= strlen(optstr); j++)
|
||||
switch (optstr[j-1])
|
||||
{
|
||||
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
|
||||
case 'd': optmesh.SplitImprove(mesh3d); break;
|
||||
case 's': optmesh.SwapImprove(mesh3d); break;
|
||||
case 'm': mesh3d.ImproveMesh(); break;
|
||||
}
|
||||
|
||||
|
||||
return MESHING3_OK;
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
MESHING3_RESULT OptimizeVolume (const MeshingParameters & mp,
|
||||
Mesh & mesh3d)
|
||||
|
||||
@ -1094,11 +1094,10 @@ static int TestSameSide (const Point3d & p1, const Point3d & p2)
|
||||
*/
|
||||
|
||||
|
||||
|
||||
void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||
void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh,
|
||||
const MeshingParameters & mp)
|
||||
{
|
||||
static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t);
|
||||
static Timer t("Mesing3::PrepareBlockFillLocalH"); RegionTimer reg(t);
|
||||
|
||||
double filldist = mp.filldist;
|
||||
|
||||
@ -1107,10 +1106,95 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||
PrintMessage (3, "blockfill local h");
|
||||
|
||||
|
||||
NgArray<Point<3> > npoints;
|
||||
|
||||
adfront -> CreateTrees();
|
||||
|
||||
double maxh = 0;
|
||||
|
||||
for (int i = 1; i <= adfront->GetNF(); i++)
|
||||
{
|
||||
const MiniElement2d & el = adfront->GetFace(i);
|
||||
for (int j = 1; j <= 3; j++)
|
||||
{
|
||||
const Point3d & p1 = adfront->GetPoint (el.PNumMod(j));
|
||||
const Point3d & p2 = adfront->GetPoint (el.PNumMod(j+1));
|
||||
|
||||
double hi = Dist (p1, p2);
|
||||
if (hi > maxh) maxh = hi;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (mp.maxh < maxh) maxh = mp.maxh;
|
||||
|
||||
// auto loch_ptr = mesh.LocalHFunction().Copy();
|
||||
// auto & loch = *loch_ptr;
|
||||
auto & loch = mesh.LocalHFunction();
|
||||
|
||||
bool changed;
|
||||
static Timer t1("loop1");
|
||||
t1.Start();
|
||||
do
|
||||
{
|
||||
loch.ClearFlags();
|
||||
|
||||
static Timer tbox("adfront-bbox");
|
||||
tbox.Start();
|
||||
for (int i = 1; i <= adfront->GetNF(); i++)
|
||||
{
|
||||
const MiniElement2d & el = adfront->GetFace(i);
|
||||
|
||||
Box<3> bbox (adfront->GetPoint (el[0]));
|
||||
bbox.Add (adfront->GetPoint (el[1]));
|
||||
bbox.Add (adfront->GetPoint (el[2]));
|
||||
|
||||
|
||||
double filld = filldist * bbox.Diam();
|
||||
bbox.Increase (filld);
|
||||
|
||||
loch.CutBoundary (bbox); // .PMin(), bbox.PMax());
|
||||
}
|
||||
tbox.Stop();
|
||||
|
||||
// locadfront = adfront;
|
||||
loch.FindInnerBoxes (adfront, NULL);
|
||||
|
||||
npoints.SetSize(0);
|
||||
loch.GetInnerPoints (npoints);
|
||||
|
||||
changed = false;
|
||||
for (int i = 1; i <= npoints.Size(); i++)
|
||||
{
|
||||
if (loch.GetH(npoints.Get(i)) > 1.5 * maxh)
|
||||
{
|
||||
loch.SetH (npoints.Get(i), maxh);
|
||||
changed = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
while (changed);
|
||||
t1.Stop();
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||
const MeshingParameters & mp)
|
||||
{
|
||||
static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t);
|
||||
|
||||
if (!mesh.HasLocalHFunction())
|
||||
{
|
||||
mesh.CalcLocalH(mp.grading);
|
||||
PrepareBlockFillLocalH(mesh, mp);
|
||||
}
|
||||
|
||||
double filldist = mp.filldist;
|
||||
|
||||
// (*testout) << "blockfill local h" << endl;
|
||||
// (*testout) << "rel filldist = " << filldist << endl;
|
||||
PrintMessage (3, "blockfill local h");
|
||||
|
||||
Box<3> bbox ( Box<3>::EMPTY_BOX );
|
||||
double maxh = 0;
|
||||
|
||||
@ -1144,44 +1228,6 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||
|
||||
if (mp.maxh < maxh) maxh = mp.maxh;
|
||||
|
||||
bool changed;
|
||||
do
|
||||
{
|
||||
mesh.LocalHFunction().ClearFlags();
|
||||
|
||||
for (int i = 1; i <= adfront->GetNF(); i++)
|
||||
{
|
||||
const MiniElement2d & el = adfront->GetFace(i);
|
||||
|
||||
Box<3> bbox (adfront->GetPoint (el[0]));
|
||||
bbox.Add (adfront->GetPoint (el[1]));
|
||||
bbox.Add (adfront->GetPoint (el[2]));
|
||||
|
||||
|
||||
double filld = filldist * bbox.Diam();
|
||||
bbox.Increase (filld);
|
||||
|
||||
mesh.LocalHFunction().CutBoundary (bbox); // .PMin(), bbox.PMax());
|
||||
}
|
||||
|
||||
// locadfront = adfront;
|
||||
mesh.LocalHFunction().FindInnerBoxes (adfront, NULL);
|
||||
|
||||
npoints.SetSize(0);
|
||||
mesh.LocalHFunction().GetInnerPoints (npoints);
|
||||
|
||||
changed = false;
|
||||
for (int i = 1; i <= npoints.Size(); i++)
|
||||
{
|
||||
if (mesh.LocalHFunction().GetH(npoints.Get(i)) > 1.5 * maxh)
|
||||
{
|
||||
mesh.LocalHFunction().SetH (npoints.Get(i), maxh);
|
||||
changed = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
while (changed);
|
||||
|
||||
if (debugparam.slowchecks)
|
||||
(*testout) << "Blockfill with points: " << endl;
|
||||
for (int i = 1; i <= npoints.Size(); i++)
|
||||
@ -1208,6 +1254,8 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||
|
||||
// find outer points
|
||||
|
||||
static Timer tloch2("build loch2");
|
||||
tloch2.Start();
|
||||
loch2.ClearFlags();
|
||||
|
||||
for (int i = 1; i <= adfront->GetNF(); i++)
|
||||
@ -1245,6 +1293,7 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||
// loch2.CutBoundary (pmin, pmax);
|
||||
loch2.CutBoundary (Box<3> (pmin, pmax)); // pmin, pmax);
|
||||
}
|
||||
tloch2.Stop();
|
||||
|
||||
// locadfront = adfront;
|
||||
loch2.FindInnerBoxes (adfront, NULL);
|
||||
|
||||
@ -28,6 +28,7 @@ class Meshing3
|
||||
NgArray<char*> problems;
|
||||
/// tolerance criterion
|
||||
double tolfak;
|
||||
NgArray<Point<3> > npoints;
|
||||
public:
|
||||
///
|
||||
Meshing3 (const string & rulefilename);
|
||||
@ -63,6 +64,7 @@ public:
|
||||
///
|
||||
void BlockFill (Mesh & mesh, double gh);
|
||||
///
|
||||
void PrepareBlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);
|
||||
void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);
|
||||
|
||||
/// uses points of adfront, and puts new elements into mesh
|
||||
|
||||
@ -852,13 +852,24 @@ namespace netgen
|
||||
{ _np = np; _typ = typ; _curved = is_curved; }
|
||||
// ar & _np & _typ & index & _curved;
|
||||
ar.DoPacked (_np, _typ, index, _curved);
|
||||
if (ar.Input())
|
||||
{ np = _np; typ = ELEMENT_TYPE(_typ); is_curved = _curved; }
|
||||
|
||||
/*
|
||||
for (size_t i = 0; i < np; i++)
|
||||
ar & pnum[i];
|
||||
*/
|
||||
if (ar.Input())
|
||||
{
|
||||
np = _np;
|
||||
typ = ELEMENT_TYPE(_typ);
|
||||
is_curved = _curved;
|
||||
flags.marked = 1;
|
||||
flags.badel = 0;
|
||||
flags.reverse = 0;
|
||||
flags.illegal = 0;
|
||||
flags.illegal_valid = 0;
|
||||
flags.badness_valid = 0;
|
||||
flags.refflag = 1;
|
||||
flags.strongrefflag = false;
|
||||
flags.deleted = 0;
|
||||
flags.fixed = 0;
|
||||
}
|
||||
|
||||
static_assert(sizeof(int) == sizeof (PointIndex));
|
||||
ar.Do( (int*)&pnum[0], np);
|
||||
}
|
||||
|
||||
@ -134,6 +134,7 @@ void ResetStatus()
|
||||
|
||||
void PushStatus(const MyStr& s)
|
||||
{
|
||||
return;
|
||||
msgstatus_stack.Append(new MyStr (s));
|
||||
SetStatMsg(s);
|
||||
threadpercent_stack.Append(0);
|
||||
@ -141,6 +142,7 @@ void PushStatus(const MyStr& s)
|
||||
|
||||
void PushStatusF(const MyStr& s)
|
||||
{
|
||||
return;
|
||||
msgstatus_stack.Append(new MyStr (s));
|
||||
SetStatMsg(s);
|
||||
threadpercent_stack.Append(0);
|
||||
@ -149,6 +151,7 @@ void PushStatusF(const MyStr& s)
|
||||
|
||||
void PopStatus()
|
||||
{
|
||||
return;
|
||||
if (msgstatus_stack.Size())
|
||||
{
|
||||
if (msgstatus_stack.Size() > 1)
|
||||
|
||||
@ -15,7 +15,6 @@ def readData(a, files):
|
||||
file=[]
|
||||
for f in files:
|
||||
for t in a[f]:
|
||||
file.append(f)
|
||||
if t['ne1d']>0:
|
||||
ne1d.append(t['ne1d'])
|
||||
if t['ne2d']>0:
|
||||
@ -24,6 +23,7 @@ def readData(a, files):
|
||||
ne3d.append(t['ne3d'])
|
||||
if t['total_badness']>0.0:
|
||||
bad.append(t['total_badness'])
|
||||
file.append(f)
|
||||
if 'angles_tet' in t:
|
||||
amin.append(t['angles_tet'][0])
|
||||
amax.append(t['angles_tet'][1])
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -87,9 +87,13 @@ def test_pickle_geom2d():
|
||||
|
||||
def test_pickle_mesh():
|
||||
import netgen.csg as csg
|
||||
geo = csg.CSGeometry()
|
||||
geo1 = csg.CSGeometry()
|
||||
geo2 = csg.CSGeometry()
|
||||
brick = csg.OrthoBrick(csg.Pnt(-3,-3,-3), csg.Pnt(3,3,3))
|
||||
mesh = geo.GenerateMesh(maxh=0.2)
|
||||
geo2.Add(brick)
|
||||
|
||||
for geo in [geo1, geo2]:
|
||||
mesh = geo.GenerateMesh(maxh=2)
|
||||
assert geo == mesh.GetGeometry()
|
||||
dump = pickle.dumps([geo,mesh])
|
||||
geo2, mesh2 = pickle.loads(dump)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user