mirror of
https://github.com/NGSolve/netgen.git
synced 2025-05-11 05:00:49 +05:00
Compare commits
No commits in common. "master" and "v6.2.2501" have entirely different histories.
@ -57,7 +57,7 @@ set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_C_COMPILER)
|
||||
set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_CXX_COMPILER)
|
||||
set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_BUILD_TYPE)
|
||||
|
||||
set(SUBPROJECT_CMAKE_ARGS "${SUBPROJECT_CMAKE_ARGS};-DCMAKE_POSITION_INDEPENDENT_CODE=ON;-DCMAKE_POLICY_VERSION_MINIMUM=3.5" CACHE INTERNAL "")
|
||||
set(SUBPROJECT_CMAKE_ARGS "${SUBPROJECT_CMAKE_ARGS};-DCMAKE_POSITION_INDEPENDENT_CODE=ON" CACHE INTERNAL "")
|
||||
|
||||
if(USE_CCACHE)
|
||||
find_program(CCACHE_FOUND NAMES ccache ccache.bat)
|
||||
|
||||
@ -109,7 +109,6 @@ if(APPLE)
|
||||
-DCMAKE_INSTALL_PREFIX=${CMAKE_INSTALL_PREFIX}/Contents/MacOS
|
||||
-DTCL_INCLUDE_PATH=${CMAKE_INSTALL_PREFIX}/Contents/Frameworks/Tcl.framework/Headers
|
||||
-DTK_INCLUDE_PATH=${CMAKE_INSTALL_PREFIX}/Contents/Frameworks/Tk.framework/Headers
|
||||
-DCMAKE_POLICY_VERSION_MINIMUM=3.5
|
||||
${SUBPROJECT_ARGS}
|
||||
)
|
||||
|
||||
|
||||
@ -718,7 +718,7 @@ namespace netgen
|
||||
|
||||
mesh -> LoadLocalMeshSize (mparam.meshsizefilename);
|
||||
for (auto mspnt : mparam.meshsize_points)
|
||||
mesh -> RestrictLocalH (mspnt.pnt, mspnt.h, mspnt.layer);
|
||||
mesh -> RestrictLocalH (mspnt.pnt, mspnt.h);
|
||||
}
|
||||
|
||||
spoints.SetSize(0);
|
||||
|
||||
@ -72,17 +72,13 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (size_t nr) const
|
||||
ret.facets.base = POINTINDEX_BASE;
|
||||
ret.facets.ptr = (int*)&el.pnum;
|
||||
|
||||
/*
|
||||
if (mesh->GetDimension() == 1)
|
||||
ret.mat = *(mesh->GetBCNamePtr(el.index-1));
|
||||
else if (mesh->GetDimension() == 2)
|
||||
ret.mat = *(mesh->GetCD2NamePtr(el.index-1));
|
||||
else
|
||||
ret.mat = *(mesh->GetCD3NamePtr(el.index-1));
|
||||
*/
|
||||
ret.mat = mesh->GetRegionName(0, el.index);
|
||||
|
||||
ret.is_curved = false;
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
@ -100,8 +96,6 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const
|
||||
ret.index = el.edgenr;
|
||||
else
|
||||
ret.index = el.si;
|
||||
|
||||
/*
|
||||
if (mesh->GetDimension() == 2)
|
||||
ret.mat = *(mesh->GetBCNamePtr(el.si-1));
|
||||
else
|
||||
@ -111,8 +105,6 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const
|
||||
else
|
||||
ret.mat = *(mesh->GetMaterialPtr(el.si));
|
||||
}
|
||||
*/
|
||||
ret.mat = mesh->GetRegionName(1, ret.index);
|
||||
|
||||
ret.points.num = el.GetNP();
|
||||
ret.points.ptr = (int*)&(el[0]);
|
||||
|
||||
@ -651,14 +651,14 @@ int Ng_FindElementOfPoint (double * p, double * lami, int build_searchtree,
|
||||
{
|
||||
Point3d p3d(p[0], p[1], p[2]);
|
||||
ind =
|
||||
mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0) + 1;
|
||||
mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0);
|
||||
}
|
||||
else
|
||||
{
|
||||
double lam3[3];
|
||||
Point3d p2d(p[0], p[1], 0);
|
||||
ind =
|
||||
mesh->GetSurfaceElementOfPoint(p2d, lam3, dummy, build_searchtree != 0) + 1;
|
||||
mesh->GetElementOfPoint(p2d, lam3, dummy, build_searchtree != 0);
|
||||
|
||||
if (ind > 0)
|
||||
{
|
||||
@ -697,7 +697,7 @@ int Ng_FindSurfaceElementOfPoint (double * p, double * lami, int build_searchtre
|
||||
{
|
||||
Point3d p3d(p[0], p[1], p[2]);
|
||||
ind =
|
||||
mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0) + 1;
|
||||
mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
@ -1011,28 +1011,68 @@ namespace netgen
|
||||
int * const indices, int numind) const
|
||||
|
||||
{
|
||||
Point<3> p(hp[0], 0., 0.);
|
||||
if(mesh->GetDimension() > 1)
|
||||
p[1] = hp[1];
|
||||
if(mesh->GetDimension() == 3)
|
||||
p[2] = hp[2];
|
||||
|
||||
for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
|
||||
switch (mesh->GetDimension())
|
||||
{
|
||||
auto & seg = (*mesh)[si];
|
||||
Point<3> p1 = (*mesh)[seg[0]];
|
||||
Point<3> p2 = (*mesh)[seg[1]];
|
||||
Vec<3> v1 = p2-p1;
|
||||
Vec<3> v2 = p-p1;
|
||||
double lam = v1*v2 / v1.Length2();
|
||||
double lam2 = (v2 - lam * v1).Length() / v1.Length();
|
||||
|
||||
if (lam >= -1e-10 && lam <= 1+1e-10 && lam2 < 1e-10)
|
||||
{
|
||||
lami[0] = 1-lam;
|
||||
return si;
|
||||
}
|
||||
case 1:
|
||||
{
|
||||
Point<3> p(hp[0], 0,0);
|
||||
for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
|
||||
{
|
||||
auto & seg = (*mesh)[si];
|
||||
Point<3> p1 = (*mesh)[seg[0]];
|
||||
Point<3> p2 = (*mesh)[seg[1]];
|
||||
double lam = (p(0)-p1(0)) / (p2(0)-p1(0));
|
||||
if (lam >= -1e-10 && lam <= 1+1e-10)
|
||||
{
|
||||
lami[0] = 1-lam;
|
||||
return si;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
case 2:
|
||||
{
|
||||
Point<3> p(hp[0], hp[1],0);
|
||||
try
|
||||
{
|
||||
auto ind = mesh->GetSurfaceElementOfPoint(p, lami, nullptr,
|
||||
build_searchtree);
|
||||
return ind - 1;
|
||||
}
|
||||
catch(const NgException & e) // quads not implemented curved yet
|
||||
{
|
||||
for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
|
||||
{
|
||||
auto & seg = (*mesh)[si];
|
||||
Point<3> p1 = (*mesh)[seg[0]];
|
||||
Point<3> p2 = (*mesh)[seg[1]];
|
||||
double lam;
|
||||
double r;
|
||||
if (fabs(p2[0]-p1[0]) >= fabs(p2[1]-p1[1]))
|
||||
{
|
||||
lam = (p[0]-p1[0])/(p2[0]-p1[0]);
|
||||
r = p[1] - p1[1] - lam*(p2[1]-p1[1]);
|
||||
}
|
||||
else
|
||||
{
|
||||
lam = (p[1]-p1[1])/(p2[1]-p1[1]);
|
||||
r = p[0] - p1[0] - lam*(p2[0]-p1[0]);
|
||||
}
|
||||
if ( lam >= -1e-10 && lam <= 1+1e-10 && fabs(r) <= 1e-10 )
|
||||
{
|
||||
lami[0] = 1-lam;
|
||||
return si;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
case 3:
|
||||
default:
|
||||
throw Exception("FindElementOfPoint<1> only implemented for mesh-dimension 1 and 2!");
|
||||
break;
|
||||
}
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
@ -1043,26 +1083,37 @@ namespace netgen
|
||||
int * const indices, int numind) const
|
||||
|
||||
{
|
||||
Point<3> pp(p[0], p[1], 0.);
|
||||
if(mesh->GetDimension() == 3)
|
||||
pp[2] = p[2];
|
||||
FlatArray<int> ind(numind, indices);
|
||||
NgArray<int> dummy(numind);
|
||||
for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
|
||||
|
||||
double lam3[3];
|
||||
auto elnr = mesh->GetSurfaceElementOfPoint(pp, lam3, ind, build_searchtree);
|
||||
if(elnr.IsValid())
|
||||
int ind;
|
||||
|
||||
if (mesh->GetDimension() == 2)
|
||||
{
|
||||
if((*mesh)[elnr].GetType() == QUAD || (*mesh)[elnr].GetType() == TRIG6)
|
||||
Point<3> p2d(p[0], p[1], 0);
|
||||
ind = mesh->GetElementOfPoint(p2d, lam3, &dummy, build_searchtree);
|
||||
}
|
||||
else
|
||||
{
|
||||
Point3d p3d(p[0], p[1], p[2]);
|
||||
ind = mesh->GetSurfaceElementOfPoint(p3d, lam3, &dummy, build_searchtree);
|
||||
}
|
||||
|
||||
if (ind > 0)
|
||||
{
|
||||
if(mesh->SurfaceElement(ind).GetType()==QUAD || mesh->SurfaceElement(ind).GetType()==TRIG6)
|
||||
{
|
||||
lami[0] = lam3[0];
|
||||
lami[1] = lam3[1];
|
||||
}
|
||||
else
|
||||
else
|
||||
{
|
||||
lami[0] = 1-lam3[0]-lam3[1];
|
||||
lami[1] = lam3[0];
|
||||
}
|
||||
}
|
||||
return elnr;
|
||||
return ind-1;
|
||||
}
|
||||
|
||||
|
||||
@ -1073,9 +1124,13 @@ namespace netgen
|
||||
int * const indices, int numind) const
|
||||
|
||||
{
|
||||
Point<3> pp(p[0], p[1], p[2]);
|
||||
FlatArray<int> ind(numind, indices);
|
||||
return mesh->GetElementOfPoint(pp, lami, ind, build_searchtree);
|
||||
NgArray<int> dummy(numind);
|
||||
for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
|
||||
|
||||
Point<3> p3d(p[0], p[1], p[2]);
|
||||
int ind =
|
||||
mesh->GetElementOfPoint(p3d, lami, &dummy, build_searchtree);
|
||||
return ind-1;
|
||||
}
|
||||
|
||||
void Ngx_Mesh :: Curve (int order)
|
||||
|
||||
@ -607,7 +607,7 @@ namespace netgen
|
||||
const_cast<Mesh&> (mesh).Compress();
|
||||
const_cast<Mesh&> (mesh).CalcSurfacesOfNode();
|
||||
const_cast<Mesh&> (mesh).RebuildSurfaceElementLists();
|
||||
const_cast<Mesh&> (mesh).BuildElementSearchTree(3);
|
||||
const_cast<Mesh&> (mesh).BuildElementSearchTree();
|
||||
|
||||
|
||||
int np = mesh.GetNP();
|
||||
|
||||
@ -66,7 +66,7 @@ void WriteFluentFormat (const Mesh & mesh,
|
||||
int i2, j2;
|
||||
NgArray<INDEX_3> surfaceelp;
|
||||
NgArray<int> surfaceeli;
|
||||
Array<ElementIndex> locels;
|
||||
NgArray<int> locels;
|
||||
|
||||
//no cells=no tets
|
||||
//no faces=2*tets
|
||||
@ -79,7 +79,7 @@ void WriteFluentFormat (const Mesh & mesh,
|
||||
snprintf (str, size(str), "(13 (4 1 %x 2 3)(",noverbface); //hexadecimal!!!
|
||||
outfile << str << endl;
|
||||
|
||||
const_cast<Mesh&> (mesh).BuildElementSearchTree(3);
|
||||
const_cast<Mesh&> (mesh).BuildElementSearchTree();
|
||||
|
||||
for (i = 1; i <= ne; i++)
|
||||
{
|
||||
@ -117,9 +117,12 @@ void WriteFluentFormat (const Mesh & mesh,
|
||||
int eli2 = 0;
|
||||
int stopsig = 0;
|
||||
|
||||
for (auto locind : locels)
|
||||
for (i2 = 1; i2 <= nel; i2++)
|
||||
{
|
||||
Element el2 = mesh[locind];
|
||||
locind = locels.Get(i2);
|
||||
//cout << " locind=" << locind << endl;
|
||||
|
||||
Element el2 = mesh.VolumeElement(locind);
|
||||
//if (inverttets)
|
||||
// el2.Invert();
|
||||
|
||||
@ -127,7 +130,7 @@ void WriteFluentFormat (const Mesh & mesh,
|
||||
{
|
||||
el2.GetFace(j2, face2);
|
||||
|
||||
if (face2.HasFace(face)) {eli2 = locind+1; stopsig = 1; break;}
|
||||
if (face2.HasFace(face)) {eli2 = locind; stopsig = 1; break;}
|
||||
}
|
||||
if (stopsig) break;
|
||||
}
|
||||
|
||||
@ -474,7 +474,7 @@ namespace netgen
|
||||
}
|
||||
|
||||
for(const auto& mspnt : mparam.meshsize_points)
|
||||
mesh.RestrictLocalH(mspnt.pnt, mspnt.h, mspnt.layer);
|
||||
mesh.RestrictLocalH(mspnt.pnt, mspnt.h);
|
||||
|
||||
mesh.LoadLocalMeshSize(mparam.meshsizefilename);
|
||||
}
|
||||
@ -1220,7 +1220,6 @@ namespace netgen
|
||||
{
|
||||
PrintMessage(3, "Optimization step ", i);
|
||||
meshopt.SetFaceIndex(k+1);
|
||||
meshopt.SetMetricWeight (mparam.elsizeweight);
|
||||
int innerstep = 0;
|
||||
for(auto optstep : mparam.optimize2d)
|
||||
{
|
||||
@ -1315,13 +1314,6 @@ namespace netgen
|
||||
if(multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES)
|
||||
return 0;
|
||||
|
||||
if(dimension == 1)
|
||||
{
|
||||
FinalizeMesh(*mesh);
|
||||
mesh->SetDimension(1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
|
||||
{
|
||||
MeshSurface(*mesh, mparam);
|
||||
|
||||
@ -4061,8 +4061,6 @@ namespace netgen
|
||||
{
|
||||
PointIndices<2> oi2(identmap[i2[0]],
|
||||
identmap[i2[1]]);
|
||||
if((!oi2[0].IsValid()) || (!oi2[1].IsValid()))
|
||||
continue;
|
||||
oi2.Sort();
|
||||
if (cutedges.Used (oi2))
|
||||
{
|
||||
|
||||
@ -14,7 +14,7 @@ namespace netgen
|
||||
|
||||
struct SpecialPointException : public Exception
|
||||
{
|
||||
SpecialPointException ()
|
||||
SpecialPointException()
|
||||
: Exception("") {}
|
||||
};
|
||||
|
||||
@ -101,14 +101,14 @@ Vec<3> CalcGrowthVector (FlatArray<Vec<3>> ns)
|
||||
return gw;
|
||||
}
|
||||
|
||||
SpecialBoundaryPoint ::GrowthGroup ::GrowthGroup (FlatArray<int> faces_,
|
||||
FlatArray<Vec<3>> normals)
|
||||
SpecialBoundaryPoint ::GrowthGroup ::GrowthGroup(FlatArray<int> faces_,
|
||||
FlatArray<Vec<3>> normals)
|
||||
{
|
||||
faces = faces_;
|
||||
growth_vector = CalcGrowthVector(normals);
|
||||
}
|
||||
|
||||
SpecialBoundaryPoint ::SpecialBoundaryPoint (
|
||||
SpecialBoundaryPoint ::SpecialBoundaryPoint(
|
||||
const std::map<int, Vec<3>>& normals)
|
||||
{
|
||||
// find opposing face normals
|
||||
@ -146,7 +146,7 @@ SpecialBoundaryPoint ::SpecialBoundaryPoint (
|
||||
growth_groups.Append(GrowthGroup(g2_faces, normals2));
|
||||
}
|
||||
|
||||
Vec<3> BoundaryLayerTool ::getEdgeTangent (PointIndex pi, int edgenr, FlatArray<Segment*> segs)
|
||||
Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr, FlatArray<Segment*> segs)
|
||||
{
|
||||
Vec<3> tangent = 0.0;
|
||||
ArrayMem<PointIndex, 2> pts;
|
||||
@ -171,7 +171,7 @@ Vec<3> BoundaryLayerTool ::getEdgeTangent (PointIndex pi, int edgenr, FlatArray<
|
||||
return tangent.Normalize();
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::LimitGrowthVectorLengths ()
|
||||
void BoundaryLayerTool ::LimitGrowthVectorLengths()
|
||||
{
|
||||
static Timer tall("BoundaryLayerTool::LimitGrowthVectorLengths");
|
||||
RegionTimer rtall(tall);
|
||||
@ -220,7 +220,7 @@ bool HaveSingleSegments (const Mesh& mesh)
|
||||
|
||||
// duplicates segments (and sets seg.si accordingly) to have a unified data
|
||||
// structure for all geometry types
|
||||
void BuildSegments (Mesh& mesh, bool have_single_segments, Array<Segment>& segments, Array<Segment>& free_segments)
|
||||
void BuildSegments (Mesh& mesh, bool have_single_segments, Array<Segment> & segments, Array<Segment> & free_segments)
|
||||
{
|
||||
// auto& topo = mesh.GetTopology();
|
||||
|
||||
@ -234,11 +234,11 @@ void BuildSegments (Mesh& mesh, bool have_single_segments, Array<Segment>& segme
|
||||
free_segments.Append(seg);
|
||||
continue;
|
||||
}
|
||||
if (!have_single_segments)
|
||||
{
|
||||
segments.Append(seg);
|
||||
continue;
|
||||
}
|
||||
if(!have_single_segments)
|
||||
{
|
||||
segments.Append(seg);
|
||||
continue;
|
||||
}
|
||||
mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els);
|
||||
for (auto seli : surf_els)
|
||||
{
|
||||
@ -284,8 +284,8 @@ void MergeAndAddSegments (Mesh& mesh, FlatArray<Segment> segments, FlatArray<Seg
|
||||
addSegment(seg);
|
||||
}
|
||||
|
||||
BoundaryLayerTool::BoundaryLayerTool (Mesh& mesh_,
|
||||
const BoundaryLayerParameters& params_)
|
||||
BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
|
||||
const BoundaryLayerParameters& params_)
|
||||
: mesh(mesh_), topo(mesh_.GetTopology()), params(params_)
|
||||
{
|
||||
static Timer timer("BoundaryLayerTool::ctor");
|
||||
@ -318,7 +318,7 @@ BoundaryLayerTool::BoundaryLayerTool (Mesh& mesh_,
|
||||
si_map[i] = i;
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::CreateNewFaceDescriptors ()
|
||||
void BoundaryLayerTool ::CreateNewFaceDescriptors()
|
||||
{
|
||||
surfacefacs.SetSize(nfd_old + 1);
|
||||
surfacefacs = 0.0;
|
||||
@ -360,7 +360,7 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors ()
|
||||
throw Exception("Surface " + to_string(si) + " is not a boundary of the domain to be grown into!");
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::CreateFaceDescriptorsSides ()
|
||||
void BoundaryLayerTool ::CreateFaceDescriptorsSides()
|
||||
{
|
||||
if (insert_only_volume_elements)
|
||||
return;
|
||||
@ -400,7 +400,7 @@ void BoundaryLayerTool ::CreateFaceDescriptorsSides ()
|
||||
}
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::CalculateGrowthVectors ()
|
||||
void BoundaryLayerTool ::CalculateGrowthVectors()
|
||||
{
|
||||
growthvectors.SetSize(np);
|
||||
growthvectors = 0.;
|
||||
@ -457,7 +457,7 @@ void BoundaryLayerTool ::CalculateGrowthVectors ()
|
||||
}
|
||||
|
||||
Array<Array<pair<SegmentIndex, int>>, SegmentIndex>
|
||||
BoundaryLayerTool ::BuildSegMap ()
|
||||
BoundaryLayerTool ::BuildSegMap()
|
||||
{
|
||||
// Bit array to keep track of segments already processed
|
||||
BitArray segs_done(nseg + 1);
|
||||
@ -536,7 +536,7 @@ BoundaryLayerTool ::BuildSegMap ()
|
||||
return segmap;
|
||||
}
|
||||
|
||||
BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface ()
|
||||
BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface()
|
||||
{
|
||||
BitArray in_surface_direction(nfd_old + 1);
|
||||
in_surface_direction.Clear();
|
||||
@ -596,7 +596,7 @@ BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface ()
|
||||
return in_surface_direction;
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::InsertNewElements (
|
||||
void BoundaryLayerTool ::InsertNewElements(
|
||||
FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap,
|
||||
const BitArray& in_surface_direction)
|
||||
{
|
||||
@ -879,61 +879,61 @@ void BoundaryLayerTool ::InsertNewElements (
|
||||
auto p2el = mesh.CreatePoint2ElementTable();
|
||||
for (SurfaceElementIndex si = 0; si < nse; si++)
|
||||
{
|
||||
// copy because surfaceels array will be resized!
|
||||
const auto sel = mesh[si];
|
||||
const auto iface = sel.GetIndex();
|
||||
|
||||
if (moved_surfaces.Test(iface))
|
||||
if (moved_surfaces.Test(sel.GetIndex()))
|
||||
{
|
||||
const auto np = sel.GetNP();
|
||||
ArrayMem<PointIndex, 4> points(sel.PNums());
|
||||
if (surfacefacs[iface] > 0)
|
||||
Array<PointIndex> points(sel.PNums());
|
||||
if (surfacefacs[sel.GetIndex()] > 0)
|
||||
Swap(points[0], points[2]);
|
||||
ArrayMem<int, 4> groups(points.Size());
|
||||
for (auto i : Range(points))
|
||||
groups[i] = getClosestGroup(points[i], si);
|
||||
groups[i] = getClosestGroup(sel[i], si);
|
||||
bool add_volume_element = true;
|
||||
for (auto pi : points)
|
||||
for (auto pi : sel.PNums())
|
||||
if (numGroups(pi) > 1)
|
||||
add_volume_element = false;
|
||||
|
||||
Element el(2 * np);
|
||||
el.PNums().Range(np, 2 * np) = points;
|
||||
auto new_index = new_mat_nrs[iface];
|
||||
if (new_index == -1)
|
||||
throw Exception("Boundary " + ToString(iface) + " with name " + mesh.GetBCName(iface - 1) + " extruded, but no new material specified for it!");
|
||||
el.SetIndex(new_index);
|
||||
|
||||
for (auto j : Range(par_heights))
|
||||
{
|
||||
el.PNums().Range(0, np) = el.PNums().Range(np, 2 * np);
|
||||
for (auto i : Range(np))
|
||||
el[np + i] = newPoint(points[i], j, groups[i]);
|
||||
auto eltype = points.Size() == 3 ? PRISM : HEX;
|
||||
Element el(eltype);
|
||||
for (auto i : Range(points))
|
||||
el[i] = points[i];
|
||||
for (auto i : Range(points))
|
||||
points[i] = newPoint(sel.PNums()[i], j, groups[i]);
|
||||
if (surfacefacs[sel.GetIndex()] > 0)
|
||||
Swap(points[0], points[2]);
|
||||
for (auto i : Range(points))
|
||||
el[sel.PNums().Size() + i] = points[i];
|
||||
auto new_index = new_mat_nrs[sel.GetIndex()];
|
||||
if (new_index == -1)
|
||||
throw Exception("Boundary " + ToString(sel.GetIndex()) + " with name " + mesh.GetBCName(sel.GetIndex() - 1) + " extruded, but no new material specified for it!");
|
||||
el.SetIndex(new_mat_nrs[sel.GetIndex()]);
|
||||
if (add_volume_element)
|
||||
mesh.AddVolumeElement(el);
|
||||
else
|
||||
{
|
||||
// Let the volume mesher fill the hole with pyramids/tets
|
||||
// To insert pyramids, we need close surface identifications on open quads
|
||||
for (auto i : Range(np))
|
||||
if (numGroups(el[i]) == 1)
|
||||
// To insert pyramids, we need close surface identifications on open
|
||||
// quads
|
||||
for (auto i : Range(points))
|
||||
if (numGroups(sel[i]) == 1)
|
||||
{
|
||||
auto pi0 = el[i];
|
||||
auto pi1 = el[np + i];
|
||||
auto pi1 = el[i + points.Size()];
|
||||
auto nr = identifications.Get(pi0, pi1);
|
||||
if (nr == 0)
|
||||
identifications.Add(pi0, pi1, identnr);
|
||||
identifications.Add(el[i], el[i + points.Size()], identnr);
|
||||
}
|
||||
}
|
||||
}
|
||||
Element2d newel = sel;
|
||||
for (auto i : Range(np))
|
||||
newel[i] = newPoint(points[i], -1, groups[i]);
|
||||
if (surfacefacs[iface] > 0)
|
||||
Swap(newel[0], newel[2]); // swap back
|
||||
newel.SetIndex(si_map[iface]);
|
||||
for (auto i : Range(points))
|
||||
newel[i] = newPoint(sel[i], -1, groups[i]);
|
||||
newel.SetIndex(si_map[sel.GetIndex()]);
|
||||
new_sels.Append(newel);
|
||||
}
|
||||
if (is_boundary_moved.Test(iface))
|
||||
if (is_boundary_moved.Test(sel.GetIndex()))
|
||||
{
|
||||
auto& sel = mesh[si];
|
||||
for (auto& p : sel.PNums())
|
||||
@ -1065,7 +1065,7 @@ void BoundaryLayerTool ::InsertNewElements (
|
||||
}
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::SetDomInOut ()
|
||||
void BoundaryLayerTool ::SetDomInOut()
|
||||
{
|
||||
if (insert_only_volume_elements)
|
||||
return;
|
||||
@ -1081,7 +1081,7 @@ void BoundaryLayerTool ::SetDomInOut ()
|
||||
}
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::SetDomInOutSides ()
|
||||
void BoundaryLayerTool ::SetDomInOutSides()
|
||||
{
|
||||
// Set the domin/domout entries for face descriptors on the "side" of new boundary layers
|
||||
if (insert_only_volume_elements)
|
||||
@ -1155,7 +1155,7 @@ void BoundaryLayerTool ::SetDomInOutSides ()
|
||||
}
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::AddSegments ()
|
||||
void BoundaryLayerTool ::AddSegments()
|
||||
{
|
||||
auto& new_segs =
|
||||
insert_only_volume_elements ? new_segments_on_moved_bnd : new_segments;
|
||||
@ -1199,14 +1199,14 @@ void BoundaryLayerTool ::AddSegments ()
|
||||
mesh.AddSegment(seg);
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::AddSurfaceElements ()
|
||||
void BoundaryLayerTool ::AddSurfaceElements()
|
||||
{
|
||||
for (auto& sel :
|
||||
insert_only_volume_elements ? new_sels_on_moved_bnd : new_sels)
|
||||
mesh.AddSurfaceElement(sel);
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::ProcessParameters ()
|
||||
void BoundaryLayerTool ::ProcessParameters()
|
||||
{
|
||||
if (int* bc = get_if<int>(¶ms.boundary); bc)
|
||||
{
|
||||
@ -1374,7 +1374,7 @@ void BoundaryLayerTool ::ProcessParameters ()
|
||||
domains.Invert();
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::Perform ()
|
||||
void BoundaryLayerTool ::Perform()
|
||||
{
|
||||
if (domains.NumSet() == 0)
|
||||
return;
|
||||
|
||||
@ -22,15 +22,15 @@ struct SpecialBoundaryPoint
|
||||
Vec<3> growth_vector;
|
||||
Array<PointIndex> new_points;
|
||||
|
||||
GrowthGroup (FlatArray<int> faces_, FlatArray<Vec<3>> normals);
|
||||
GrowthGroup (const GrowthGroup&) = default;
|
||||
GrowthGroup () = default;
|
||||
GrowthGroup(FlatArray<int> faces_, FlatArray<Vec<3>> normals);
|
||||
GrowthGroup(const GrowthGroup&) = default;
|
||||
GrowthGroup() = default;
|
||||
};
|
||||
Array<GrowthGroup> growth_groups;
|
||||
Vec<3> separating_direction;
|
||||
|
||||
SpecialBoundaryPoint (const std::map<int, Vec<3>>& normals);
|
||||
SpecialBoundaryPoint () = default;
|
||||
SpecialBoundaryPoint(const std::map<int, Vec<3>>& normals);
|
||||
SpecialBoundaryPoint() = default;
|
||||
};
|
||||
|
||||
DLL_HEADER void GenerateBoundaryLayer (Mesh& mesh,
|
||||
@ -41,7 +41,7 @@ DLL_HEADER int /* new_domain_number */ GenerateBoundaryLayer2 (Mesh& mesh, int d
|
||||
class BoundaryLayerTool
|
||||
{
|
||||
public:
|
||||
BoundaryLayerTool (Mesh& mesh_, const BoundaryLayerParameters& params_);
|
||||
BoundaryLayerTool(Mesh& mesh_, const BoundaryLayerParameters& params_);
|
||||
void ProcessParameters ();
|
||||
void Perform ();
|
||||
|
||||
|
||||
@ -243,6 +243,10 @@ namespace netgen
|
||||
|
||||
Array<SegmentIndex> segments;
|
||||
|
||||
// surface index map
|
||||
Array<int> si_map(mesh.GetNFD()+2);
|
||||
si_map = -1;
|
||||
|
||||
// int fd_old = mesh.GetNFD();
|
||||
|
||||
int max_edge_nr = -1;
|
||||
@ -250,8 +254,8 @@ namespace netgen
|
||||
|
||||
for(const auto& seg : line_segments)
|
||||
{
|
||||
if(seg.edgenr > max_edge_nr)
|
||||
max_edge_nr = seg.edgenr;
|
||||
if(seg.epgeominfo[0].edgenr > max_edge_nr)
|
||||
max_edge_nr = seg.epgeominfo[0].edgenr;
|
||||
if(seg.domin > max_domain)
|
||||
max_domain = seg.domin;
|
||||
if(seg.domout > max_domain)
|
||||
@ -259,7 +263,6 @@ namespace netgen
|
||||
}
|
||||
|
||||
int new_domain = max_domain+1;
|
||||
int new_edge_nr = max_edge_nr+1;
|
||||
|
||||
BitArray active_boundaries(max_edge_nr+1);
|
||||
BitArray active_segments(nseg);
|
||||
@ -285,10 +288,7 @@ namespace netgen
|
||||
// int new_fd_index =
|
||||
mesh.AddFaceDescriptor(new_fd);
|
||||
if(should_make_new_domain)
|
||||
{
|
||||
mesh.SetMaterial(new_domain, "layer_" + mesh.GetMaterial(domain));
|
||||
mesh.SetBCName(new_edge_nr - 1, "moved");
|
||||
}
|
||||
mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(domain-1));
|
||||
}
|
||||
|
||||
for(auto segi : Range(line_segments))
|
||||
@ -618,6 +618,8 @@ namespace netgen
|
||||
}
|
||||
}
|
||||
|
||||
map<pair<PointIndex, PointIndex>, int> seg2edge;
|
||||
|
||||
// insert new elements ( and move old ones )
|
||||
for(auto si : moved_segs)
|
||||
{
|
||||
@ -627,6 +629,8 @@ namespace netgen
|
||||
auto & pm0 = mapto[seg[0]];
|
||||
auto & pm1 = mapto[seg[1]];
|
||||
|
||||
// auto newindex = si_map[domain];
|
||||
|
||||
Segment s = seg;
|
||||
s.geominfo[0] = {};
|
||||
s.geominfo[1] = {};
|
||||
@ -634,10 +638,10 @@ namespace netgen
|
||||
s[1] = pm1.Last();
|
||||
s[2] = PointIndex::INVALID;
|
||||
auto pair = s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]);
|
||||
s.edgenr = new_edge_nr;
|
||||
s.epgeominfo[0].edgenr = -1;
|
||||
s.epgeominfo[1].edgenr = -1;
|
||||
s.si = s.edgenr;
|
||||
if(seg2edge.find(pair) == seg2edge.end())
|
||||
seg2edge[pair] = ++max_edge_nr;
|
||||
s.edgenr = seg2edge[pair];
|
||||
s.si = seg.si;
|
||||
mesh.AddSegment(s);
|
||||
|
||||
for ( auto i : Range(thicknesses))
|
||||
|
||||
@ -70,7 +70,7 @@ BuildNeighbors (FlatArray<PointIndex> points, const Mesh& mesh)
|
||||
return neighbors;
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::InterpolateGrowthVectors ()
|
||||
void BoundaryLayerTool ::InterpolateGrowthVectors()
|
||||
{
|
||||
point_types.SetSize(mesh.GetNP());
|
||||
for (auto p : mesh.Points().Range())
|
||||
@ -139,7 +139,7 @@ void BoundaryLayerTool ::InterpolateGrowthVectors ()
|
||||
faces.Append(sei);
|
||||
}
|
||||
|
||||
if (faces.Size() == 2 && mesh[faces[0]].GetIndex() != mesh[faces[1]].GetIndex())
|
||||
if (faces.Size() == 2)
|
||||
{
|
||||
auto n0 = getNormal(mesh[faces[0]]);
|
||||
auto n1 = getNormal(mesh[faces[1]]);
|
||||
@ -299,7 +299,7 @@ void BoundaryLayerTool ::InterpolateGrowthVectors ()
|
||||
}
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors ()
|
||||
void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors()
|
||||
{
|
||||
static Timer tall("InterpolateSurfaceGrowthVectors");
|
||||
RegionTimer rtall(tall);
|
||||
@ -399,7 +399,7 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors ()
|
||||
growthvectors[pi] += corrections[pi];
|
||||
}
|
||||
|
||||
void BoundaryLayerTool ::FixSurfaceElements ()
|
||||
void BoundaryLayerTool ::FixSurfaceElements()
|
||||
{
|
||||
static Timer tall("FixSurfaceElements");
|
||||
RegionTimer rtall(tall);
|
||||
|
||||
@ -29,7 +29,7 @@ struct GrowthVectorLimiter
|
||||
Array<PointIndex, PointIndex> map_from;
|
||||
Table<SurfaceElementIndex, PointIndex> p2sel;
|
||||
|
||||
GrowthVectorLimiter (BoundaryLayerTool& tool_)
|
||||
GrowthVectorLimiter(BoundaryLayerTool& tool_)
|
||||
: tool(tool_), params(tool_.params), mesh(tool_.mesh), height(tool_.total_height), growthvectors(tool_.growthvectors), map_from(mesh.Points().Size())
|
||||
{
|
||||
changed_domains = tool.domains;
|
||||
@ -287,7 +287,7 @@ struct GrowthVectorLimiter
|
||||
if (factor == 0.0)
|
||||
return;
|
||||
// for (PointIndex pi : IntRange(tool.np, mesh.GetNP()))
|
||||
for (PointIndex pi : mesh.Points().Range().Modify(tool.np, 0))
|
||||
for (PointIndex pi : mesh.Points().Range().Modify(tool.np,0))
|
||||
{
|
||||
// auto pi_from = map_from[pi];
|
||||
std::set<PointIndex> pis;
|
||||
@ -464,7 +464,8 @@ struct GrowthVectorLimiter
|
||||
}
|
||||
|
||||
template <typename TFunc>
|
||||
void FindTreeIntersections (double trig_shift, double seg_shift, TFunc f, TBitArray<PointIndex>* relevant_points = nullptr)
|
||||
void FindTreeIntersections (double trig_shift, double seg_shift, TFunc f,
|
||||
TBitArray<PointIndex> * relevant_points = nullptr)
|
||||
{
|
||||
static Timer t("GrowthVectorLimiter::FindTreeIntersections");
|
||||
RegionTimer rt(t);
|
||||
@ -505,25 +506,6 @@ struct GrowthVectorLimiter
|
||||
RegionTimer reg(t);
|
||||
// check if surface trigs are intersecting each other
|
||||
bool changed = true;
|
||||
std::set<PointIndex> special_points;
|
||||
|
||||
if (tool.insert_only_volume_elements)
|
||||
for (auto [pi, special_point] : tool.special_boundary_points)
|
||||
{
|
||||
special_points.insert(pi);
|
||||
for (auto& group : special_point.growth_groups)
|
||||
special_points.insert(group.new_points.Last());
|
||||
}
|
||||
|
||||
auto skip_trig = [&] (const Element2d& tri) {
|
||||
if (!tool.insert_only_volume_elements)
|
||||
return false;
|
||||
for (auto pi : tri.PNums())
|
||||
if (special_points.find(pi) != special_points.end())
|
||||
return true;
|
||||
return false;
|
||||
};
|
||||
|
||||
while (changed)
|
||||
{
|
||||
changed = false;
|
||||
@ -535,9 +517,6 @@ struct GrowthVectorLimiter
|
||||
{
|
||||
const Element2d& tri = Get(sei);
|
||||
|
||||
if (skip_trig(tri))
|
||||
continue;
|
||||
|
||||
Box<3> box(Box<3>::EMPTY_BOX);
|
||||
for (PointIndex pi : tri.PNums())
|
||||
box.Add(GetPoint(pi, 1.0, true));
|
||||
@ -550,9 +529,6 @@ struct GrowthVectorLimiter
|
||||
{
|
||||
const Element2d& tri = Get(sei);
|
||||
|
||||
if (skip_trig(tri))
|
||||
continue;
|
||||
|
||||
Box<3> box(Box<3>::EMPTY_BOX);
|
||||
for (PointIndex pi : tri.PNums())
|
||||
box.Add(GetPoint(pi, 1.0, true));
|
||||
@ -653,7 +629,7 @@ struct GrowthVectorLimiter
|
||||
size_t limit_counter = 1;
|
||||
|
||||
TBitArray<PointIndex> relevant_points, relevant_points_next;
|
||||
relevant_points.SetSize(mesh.Points().Size() + 1);
|
||||
relevant_points.SetSize(mesh.Points().Size() + 1);
|
||||
relevant_points_next.SetSize(mesh.Points().Size() + 1);
|
||||
relevant_points.Set();
|
||||
|
||||
@ -709,9 +685,8 @@ struct GrowthVectorLimiter
|
||||
for (auto pi : Range(growthvectors))
|
||||
check_point(pi);
|
||||
|
||||
if (!tool.insert_only_volume_elements)
|
||||
for (auto& [special_pi, special_point] : tool.special_boundary_points)
|
||||
check_point(special_pi);
|
||||
for (auto& [special_pi, special_point] : tool.special_boundary_points)
|
||||
check_point(special_pi);
|
||||
}
|
||||
|
||||
void Perform ()
|
||||
|
||||
@ -9,25 +9,6 @@
|
||||
|
||||
namespace netgen
|
||||
{
|
||||
inline int GetVolElement(const Mesh& mesh, const Point<3>& p,
|
||||
double* lami)
|
||||
{
|
||||
if(mesh.GetDimension() == 3)
|
||||
{
|
||||
auto ei = mesh.GetElementOfPoint(p, lami, true);
|
||||
if(!ei.IsValid())
|
||||
return -1;
|
||||
return ei;
|
||||
}
|
||||
else
|
||||
{
|
||||
auto ei = mesh.GetSurfaceElementOfPoint(p, lami, true);
|
||||
if(!ei.IsValid())
|
||||
return -1;
|
||||
return ei;
|
||||
}
|
||||
}
|
||||
|
||||
RKStepper :: ~RKStepper()
|
||||
{
|
||||
delete a;
|
||||
@ -209,9 +190,9 @@ namespace netgen
|
||||
|
||||
for(int i=0; i<potential_startpoints.Size(); i++)
|
||||
{
|
||||
int elnr = GetVolElement(mesh, potential_startpoints[i], lami);
|
||||
if (elnr == -1)
|
||||
continue;
|
||||
int elnr = mesh.GetElementOfPoint(potential_startpoints[i],lami,true) - 1;
|
||||
if(elnr == -1)
|
||||
continue;
|
||||
|
||||
mesh.SetPointSearchStartElement(elnr);
|
||||
|
||||
@ -308,7 +289,8 @@ namespace netgen
|
||||
dirstart.SetSize(0);
|
||||
dirstart.Append(0);
|
||||
|
||||
int startelnr = GetVolElement(mesh, startpoint,startlami);
|
||||
|
||||
int startelnr = mesh.GetElementOfPoint(startpoint,startlami,true) - 1;
|
||||
(*testout) << "p = " << startpoint << "; elnr = " << startelnr << endl;
|
||||
if (startelnr == -1)
|
||||
return;
|
||||
@ -360,7 +342,7 @@ namespace netgen
|
||||
Point<3> newp;
|
||||
while(!stepper.GetNextData(newp,dummyt,h) && elnr != -1)
|
||||
{
|
||||
elnr = GetVolElement(mesh, newp, lami);
|
||||
elnr = mesh.GetElementOfPoint(newp,lami,true) - 1;
|
||||
if(elnr != -1)
|
||||
{
|
||||
mesh.SetPointSearchStartElement(elnr);
|
||||
|
||||
@ -2194,7 +2194,7 @@ void MeshOptimize3d :: SwapImproveSurface (
|
||||
|
||||
double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
|
||||
Table<ElementIndex, PointIndex> & elementsonnode,
|
||||
DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool conform_segments, bool check_only )
|
||||
DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool check_only )
|
||||
{
|
||||
PointIndex pi1, pi2, pi3, pi4, pi5;
|
||||
Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
|
||||
@ -2228,31 +2228,28 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
|
||||
}
|
||||
|
||||
|
||||
if(!conform_segments)
|
||||
{
|
||||
bool bface = 0;
|
||||
for (int k = 0; k < belementsonnode[pi1].Size(); k++)
|
||||
{
|
||||
const Element2d & bel =
|
||||
mesh[belementsonnode[pi1][k]];
|
||||
bool bface = 0;
|
||||
for (int k = 0; k < belementsonnode[pi1].Size(); k++)
|
||||
{
|
||||
const Element2d & bel =
|
||||
mesh[belementsonnode[pi1][k]];
|
||||
|
||||
bool bface1 = 1;
|
||||
for (int l = 0; l < 3; l++)
|
||||
if (bel[l] != pi1 && bel[l] != pi2 && bel[l] != pi3)
|
||||
{
|
||||
bface1 = 0;
|
||||
break;
|
||||
}
|
||||
|
||||
if (bface1)
|
||||
bool bface1 = 1;
|
||||
for (int l = 0; l < 3; l++)
|
||||
if (bel[l] != pi1 && bel[l] != pi2 && bel[l] != pi3)
|
||||
{
|
||||
bface = 1;
|
||||
bface1 = 0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (bface) return 0.0;
|
||||
}
|
||||
if (bface1)
|
||||
{
|
||||
bface = 1;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (bface) return 0.0;
|
||||
|
||||
|
||||
FlatArray<ElementIndex> row = elementsonnode[pi1];
|
||||
@ -2370,11 +2367,11 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
|
||||
2 -> 3 conversion
|
||||
*/
|
||||
|
||||
void MeshOptimize3d :: SwapImprove2 (bool conform_segments)
|
||||
void MeshOptimize3d :: SwapImprove2 ()
|
||||
{
|
||||
static Timer t("MeshOptimize3d::SwapImprove2"); RegionTimer reg(t);
|
||||
|
||||
if (!conform_segments && goal == OPT_CONFORM) return;
|
||||
if (goal == OPT_CONFORM) return;
|
||||
|
||||
mesh.BuildBoundaryEdges(false);
|
||||
|
||||
@ -2436,7 +2433,7 @@ void MeshOptimize3d :: SwapImprove2 (bool conform_segments)
|
||||
|
||||
for (int j = 0; j < 4; j++)
|
||||
{
|
||||
double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, conform_segments, true);
|
||||
double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, true);
|
||||
if(d_badness<0.0)
|
||||
my_faces_with_improvement.Append( std::make_tuple(d_badness, eli1, j) );
|
||||
}
|
||||
@ -2452,7 +2449,7 @@ void MeshOptimize3d :: SwapImprove2 (bool conform_segments)
|
||||
{
|
||||
if(mesh[eli].IsDeleted())
|
||||
continue;
|
||||
if(SwapImprove2( eli, j, elementsonnode, belementsonnode, conform_segments, false) < 0.0)
|
||||
if(SwapImprove2( eli, j, elementsonnode, belementsonnode, false) < 0.0)
|
||||
cnt++;
|
||||
}
|
||||
|
||||
|
||||
@ -45,8 +45,8 @@ public:
|
||||
void SwapImprove (const TBitArray<ElementIndex> * working_elements = NULL);
|
||||
void SwapImproveSurface (const TBitArray<ElementIndex> * working_elements = NULL,
|
||||
const NgArray< idmap_type* > * idmaps = NULL);
|
||||
void SwapImprove2 (bool conform_segments = false);
|
||||
double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool conform_segments, bool check_only=false );
|
||||
void SwapImprove2 ();
|
||||
double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool check_only=false );
|
||||
|
||||
void ImproveMesh() { mesh.ImproveMesh(mp, goal); }
|
||||
|
||||
|
||||
@ -196,41 +196,54 @@ namespace netgen
|
||||
fabs (p(1) - root->xmid[1]) > root->h2)
|
||||
return;
|
||||
|
||||
GradingBox * box = Find(p);
|
||||
if (box->HOpt() <= 1.2 * h) return;
|
||||
if (GetH(p) <= 1.2 * h) return;
|
||||
|
||||
GradingBox * box = root;
|
||||
GradingBox * nbox = root;
|
||||
GradingBox * ngb;
|
||||
int childnr;
|
||||
double x1[3], x2[3];
|
||||
|
||||
while (nbox)
|
||||
{
|
||||
box = nbox;
|
||||
childnr = 0;
|
||||
if (p(0) > box->xmid[0]) childnr += 1;
|
||||
if (p(1) > box->xmid[1]) childnr += 2;
|
||||
nbox = box->childs[childnr];
|
||||
};
|
||||
|
||||
while (2 * box->h2 > h)
|
||||
{
|
||||
int childnr = 0;
|
||||
double x1[3], x2[3];
|
||||
|
||||
childnr = 0;
|
||||
if (p(0) > box->xmid[0]) childnr += 1;
|
||||
if (p(1) > box->xmid[1]) childnr += 2;
|
||||
|
||||
double h2 = box->h2;
|
||||
if (p(0) > box->xmid[0])
|
||||
if (childnr & 1)
|
||||
{
|
||||
childnr += 1;
|
||||
x1[0] = box->xmid[0];
|
||||
x2[0] = x1[0]+h2;
|
||||
x2[0] = x1[0]+h2; // box->x2[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
x2[0] = box->xmid[0];
|
||||
x1[0] = x2[0]-h2;
|
||||
x1[0] = x2[0]-h2; // box->x1[0];
|
||||
}
|
||||
|
||||
if (p(1) > box->xmid[1])
|
||||
if (childnr & 2)
|
||||
{
|
||||
childnr += 2;
|
||||
x1[1] = box->xmid[1];
|
||||
x2[1] = x1[1]+h2;
|
||||
x2[1] = x1[1]+h2; // box->x2[1];
|
||||
}
|
||||
else
|
||||
{
|
||||
x2[1] = box->xmid[1];
|
||||
x1[1] = x2[1]-h2;
|
||||
x1[1] = x2[1]-h2; // box->x1[1];
|
||||
}
|
||||
x1[2] = x2[2] = 0;
|
||||
|
||||
auto ngb = new GradingBox (x1, x2);
|
||||
ngb = new GradingBox (x1, x2);
|
||||
box->childs[childnr] = ngb;
|
||||
ngb->father = box;
|
||||
|
||||
@ -263,52 +276,67 @@ namespace netgen
|
||||
fabs (p(2) - root->xmid[2]) > root->h2)
|
||||
return;
|
||||
|
||||
GradingBox * box = Find(p);
|
||||
if (box->HOpt() <= 1.2 * h) return;
|
||||
if (GetH(p) <= 1.2 * h) return;
|
||||
|
||||
GradingBox * box = root;
|
||||
GradingBox * nbox = root;
|
||||
GradingBox * ngb;
|
||||
int childnr;
|
||||
double x1[3], x2[3];
|
||||
|
||||
while (nbox)
|
||||
{
|
||||
box = nbox;
|
||||
childnr = 0;
|
||||
if (p(0) > box->xmid[0]) childnr += 1;
|
||||
if (p(1) > box->xmid[1]) childnr += 2;
|
||||
if (p(2) > box->xmid[2]) childnr += 4;
|
||||
nbox = box->childs[childnr];
|
||||
};
|
||||
|
||||
|
||||
while (2 * box->h2 > h)
|
||||
{
|
||||
double x1[3], x2[3];
|
||||
int childnr = 0;
|
||||
childnr = 0;
|
||||
if (p(0) > box->xmid[0]) childnr += 1;
|
||||
if (p(1) > box->xmid[1]) childnr += 2;
|
||||
if (p(2) > box->xmid[2]) childnr += 4;
|
||||
|
||||
double h2 = box->h2;
|
||||
|
||||
if (p(0) > box->xmid[0])
|
||||
if (childnr & 1)
|
||||
{
|
||||
childnr += 1;
|
||||
x1[0] = box->xmid[0];
|
||||
x2[0] = x1[0]+h2;
|
||||
x2[0] = x1[0]+h2; // box->x2[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
x2[0] = box->xmid[0];
|
||||
x1[0] = x2[0]-h2;
|
||||
x1[0] = x2[0]-h2; // box->x1[0];
|
||||
}
|
||||
|
||||
if (p(1) > box->xmid[1])
|
||||
if (childnr & 2)
|
||||
{
|
||||
childnr += 2;
|
||||
x1[1] = box->xmid[1];
|
||||
x2[1] = x1[1]+h2;
|
||||
x2[1] = x1[1]+h2; // box->x2[1];
|
||||
}
|
||||
else
|
||||
{
|
||||
x2[1] = box->xmid[1];
|
||||
x1[1] = x2[1]-h2;
|
||||
x1[1] = x2[1]-h2; // box->x1[1];
|
||||
}
|
||||
|
||||
if (p(2) > box->xmid[2])
|
||||
if (childnr & 4)
|
||||
{
|
||||
childnr += 4;
|
||||
x1[2] = box->xmid[2];
|
||||
x2[2] = x1[2]+h2;
|
||||
x2[2] = x1[2]+h2; // box->x2[2];
|
||||
}
|
||||
else
|
||||
{
|
||||
x2[2] = box->xmid[2];
|
||||
x1[2] = x2[2]-h2;
|
||||
x1[2] = x2[2]-h2; // box->x1[2];
|
||||
}
|
||||
|
||||
auto ngb = new GradingBox (x1, x2);
|
||||
ngb = new GradingBox (x1, x2);
|
||||
box->childs[childnr] = ngb;
|
||||
ngb->father = box;
|
||||
|
||||
@ -338,7 +366,37 @@ namespace netgen
|
||||
|
||||
double LocalH :: GetH (Point<3> x) const
|
||||
{
|
||||
return Find(x)->HOpt();
|
||||
const GradingBox * box = root;
|
||||
if (dimension == 2)
|
||||
{
|
||||
while (1)
|
||||
{
|
||||
int childnr = 0;
|
||||
if (x(0) > box->xmid[0]) childnr += 1;
|
||||
if (x(1) > box->xmid[1]) childnr += 2;
|
||||
|
||||
if (box->childs[childnr])
|
||||
box = box->childs[childnr];
|
||||
else
|
||||
return box->hopt;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
while (1)
|
||||
{
|
||||
int childnr = 0;
|
||||
if (x(0) > box->xmid[0]) childnr += 1;
|
||||
if (x(1) > box->xmid[1]) childnr += 2;
|
||||
if (x(2) > box->xmid[2]) childnr += 4;
|
||||
|
||||
if (box->childs[childnr])
|
||||
box = box->childs[childnr];
|
||||
else
|
||||
return box->hopt;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
@ -430,41 +488,6 @@ namespace netgen
|
||||
}
|
||||
|
||||
|
||||
GradingBox * LocalH :: Find (Point<3> p) const
|
||||
{
|
||||
GradingBox * box = root;
|
||||
if (dimension == 2)
|
||||
{
|
||||
while (1)
|
||||
{
|
||||
int childnr = 0;
|
||||
if (p(0) > box->xmid[0]) childnr += 1;
|
||||
if (p(1) > box->xmid[1]) childnr += 2;
|
||||
|
||||
if (box->childs[childnr])
|
||||
box = box->childs[childnr];
|
||||
else
|
||||
return box;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
while (1)
|
||||
{
|
||||
int childnr = 0;
|
||||
if (p(0) > box->xmid[0]) childnr += 1;
|
||||
if (p(1) > box->xmid[1]) childnr += 2;
|
||||
if (p(2) > box->xmid[2]) childnr += 4;
|
||||
|
||||
if (box->childs[childnr])
|
||||
box = box->childs[childnr];
|
||||
else
|
||||
return box;
|
||||
}
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
|
||||
void LocalH :: FindInnerBoxes (const AdFront3 & adfront,
|
||||
int (*testinner)(const Point3d & p1))
|
||||
|
||||
@ -54,7 +54,6 @@ namespace netgen
|
||||
|
||||
Point<3> PMid() const { return Point<3> (xmid[0], xmid[1], xmid[2]); }
|
||||
double H2() const { return h2; }
|
||||
double HOpt() const { return hopt; }
|
||||
|
||||
bool HasChilds() const
|
||||
{
|
||||
@ -120,8 +119,6 @@ namespace netgen
|
||||
|
||||
void CutBoundary (const Box<3> & box)
|
||||
{ CutBoundaryRec (box.PMin(), box.PMax(), root); }
|
||||
|
||||
GradingBox * Find(Point<3> p) const;
|
||||
|
||||
/// find inner boxes
|
||||
void FindInnerBoxes (const class AdFront3 & adfront,
|
||||
|
||||
@ -12,15 +12,15 @@
|
||||
|
||||
namespace netgen
|
||||
{
|
||||
ElementIndex Find3dElement (const Mesh& mesh,
|
||||
const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
optional<FlatArray<int>> indices,
|
||||
BoxTree<3, ElementIndex> * searchtree,
|
||||
const bool allowindex)
|
||||
int Find3dElement (const Mesh& mesh,
|
||||
const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
const NgArray<int> * const indices,
|
||||
BoxTree<3> * searchtree,
|
||||
const bool allowindex = true)
|
||||
{
|
||||
int ne = 0;
|
||||
Array<ElementIndex> locels;
|
||||
NgArray<int> locels;
|
||||
if (searchtree)
|
||||
{
|
||||
searchtree->GetIntersecting (p, p, locels);
|
||||
@ -29,90 +29,92 @@ namespace netgen
|
||||
else
|
||||
ne = mesh.GetNE();
|
||||
|
||||
for (auto i : Range(ne))
|
||||
for (int i = 1; i <= ne; i++)
|
||||
{
|
||||
ElementIndex ei;
|
||||
int ii;
|
||||
|
||||
if (searchtree)
|
||||
ei = locels[i];
|
||||
ii = locels.Get(i);
|
||||
else
|
||||
ei = i;
|
||||
ii = i;
|
||||
|
||||
if(indices && indices->Size() > 0)
|
||||
if(indices != NULL && indices->Size() > 0)
|
||||
{
|
||||
bool contained = indices->Contains(mesh[ei].GetIndex());
|
||||
bool contained = indices->Contains(mesh.VolumeElement(ii).GetIndex());
|
||||
if((allowindex && !contained) || (!allowindex && contained)) continue;
|
||||
}
|
||||
|
||||
if(mesh.PointContainedIn3DElement(p,lami,ei))
|
||||
return ei;
|
||||
if(mesh.PointContainedIn3DElement(p,lami,ii))
|
||||
return ii;
|
||||
}
|
||||
|
||||
// Not found, try uncurved variant:
|
||||
for (auto i : Range(ne))
|
||||
for (int i = 1; i <= ne; i++)
|
||||
{
|
||||
ElementIndex ei;
|
||||
int ii;
|
||||
|
||||
if (searchtree)
|
||||
ei = locels[i];
|
||||
ii = locels.Get(i);
|
||||
else
|
||||
ei = i;
|
||||
ii = i;
|
||||
|
||||
if(indices && indices->Size() > 0)
|
||||
if(indices != NULL && indices->Size() > 0)
|
||||
{
|
||||
bool contained = indices->Contains(mesh[ei].GetIndex());
|
||||
bool contained = indices->Contains(mesh.VolumeElement(ii).GetIndex());
|
||||
if((allowindex && !contained) || (!allowindex && contained)) continue;
|
||||
}
|
||||
|
||||
|
||||
if(mesh.PointContainedIn3DElementOld(p,lami,ei))
|
||||
if(mesh.PointContainedIn3DElementOld(p,lami,ii))
|
||||
{
|
||||
(*testout) << "WARNING: found element of point " << p <<" only for uncurved mesh" << endl;
|
||||
return ei;
|
||||
return ii;
|
||||
}
|
||||
}
|
||||
return ElementIndex::INVALID;
|
||||
return 0;
|
||||
}
|
||||
|
||||
SurfaceElementIndex
|
||||
Find2dElement (const Mesh& mesh,
|
||||
const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
std::optional<FlatArray<int>> indices,
|
||||
BoxTree<3, SurfaceElementIndex> * searchtree,
|
||||
bool allowindex)
|
||||
int Find2dElement (const Mesh& mesh,
|
||||
const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
const NgArray<int> * const indices,
|
||||
BoxTree<3> * searchtree,
|
||||
const bool allowindex = true)
|
||||
{
|
||||
double vlam[3];
|
||||
ElementIndex velement = ElementIndex::INVALID;
|
||||
int velement = 0;
|
||||
|
||||
if(mesh.GetNE())
|
||||
{
|
||||
if(searchtree)
|
||||
const_cast<Mesh&>(mesh).BuildElementSearchTree(3);
|
||||
velement = Find3dElement(mesh, p,vlam, nullopt,searchtree ? mesh.GetElementSearchTree() : nullptr,allowindex);
|
||||
}
|
||||
velement = Find3dElement(mesh, p,vlam,NULL,searchtree,allowindex);
|
||||
|
||||
//(*testout) << "p " << p << endl;
|
||||
//(*testout) << "velement " << velement << endl;
|
||||
|
||||
// first try to find a volume element containing p and project to face
|
||||
if(velement.IsValid())
|
||||
if(velement!=0)
|
||||
{
|
||||
auto & topology = mesh.GetTopology();
|
||||
const auto & fnrs = topology.GetFaces(velement);
|
||||
auto faces = ArrayMem<SurfaceElementIndex,4>();
|
||||
for(auto face : fnrs)
|
||||
faces.Append(topology.GetFace2SurfaceElement(face));
|
||||
// NgArray<int> faces;
|
||||
// topology.GetElementFaces(velement,faces);
|
||||
auto faces = Array<int> (topology.GetFaces(ElementIndex(velement-1)));
|
||||
|
||||
//(*testout) << "faces " << faces << endl;
|
||||
|
||||
for(int i=0; i<faces.Size(); i++)
|
||||
faces[i] = topology.GetFace2SurfaceElement(faces[i])+1;
|
||||
|
||||
//(*testout) << "surfel " << faces << endl;
|
||||
|
||||
for(int i=0; i<faces.Size(); i++)
|
||||
{
|
||||
if(!faces[i].IsValid())
|
||||
if(faces[i] == 0)
|
||||
continue;
|
||||
auto sel = mesh.SurfaceElement(faces[i]);
|
||||
if(indices && indices->Size() > 0 && !indices->Contains(sel.GetIndex()))
|
||||
if(indices && indices->Size() != 0 && !indices->Contains(sel.GetIndex()))
|
||||
continue;
|
||||
|
||||
auto & el = mesh[velement];
|
||||
auto & el = mesh.VolumeElement(velement);
|
||||
|
||||
if (el.GetType() == TET)
|
||||
{
|
||||
double lam4[4] = { vlam[0], vlam[1], vlam[2], 1.0-vlam[0]-vlam[1]-vlam[2] };
|
||||
@ -125,7 +127,7 @@ namespace netgen
|
||||
for(auto k : Range(4))
|
||||
if(sel[j] == el[k])
|
||||
lami[j-1] = lam4[k]/(1.0-face_lam);
|
||||
return SurfaceElementIndex(faces[i]);
|
||||
return faces[i];
|
||||
}
|
||||
}
|
||||
|
||||
@ -137,8 +139,9 @@ namespace netgen
|
||||
// Did't find any matching face of a volume element, search 2d elements directly
|
||||
int ne;
|
||||
|
||||
Array<SurfaceElementIndex> locels;
|
||||
if (searchtree)
|
||||
NgArray<int> locels;
|
||||
// TODO: build search tree for surface elements
|
||||
if (!mesh.GetNE() && searchtree)
|
||||
{
|
||||
searchtree->GetIntersecting (p, p, locels);
|
||||
ne = locels.Size();
|
||||
@ -146,38 +149,37 @@ namespace netgen
|
||||
else
|
||||
ne = mesh.GetNSE();
|
||||
|
||||
for (auto i : Range(ne))
|
||||
for (int i = 1; i <= ne; i++)
|
||||
{
|
||||
SurfaceElementIndex ii;
|
||||
int ii;
|
||||
|
||||
if (locels.Size())
|
||||
ii = locels[i];
|
||||
ii = locels.Get(i);
|
||||
else
|
||||
ii = i;
|
||||
|
||||
if(indices && indices->Size() > 0)
|
||||
if(indices != NULL && indices->Size() > 0)
|
||||
{
|
||||
bool contained = indices->Contains(mesh[ii].GetIndex());
|
||||
bool contained = indices->Contains(mesh.SurfaceElement(ii).GetIndex());
|
||||
if((allowindex && !contained) || (!allowindex && contained)) continue;
|
||||
}
|
||||
if(mesh.PointContainedIn2DElement(p,lami,ii))
|
||||
return ii;
|
||||
|
||||
if(mesh.PointContainedIn2DElement(p,lami,ii)) return ii;
|
||||
|
||||
}
|
||||
return SurfaceElementIndex::INVALID;
|
||||
return 0;
|
||||
}
|
||||
|
||||
SegmentIndex Find1dElement (const Mesh& mesh,
|
||||
const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
std::optional<FlatArray<int>> indices,
|
||||
BoxTree<3> * searchtree,
|
||||
const bool allowindex = true)
|
||||
int Find1dElement (const Mesh& mesh,
|
||||
const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
const NgArray<int> * const indices,
|
||||
BoxTree<3> * searchtree,
|
||||
const bool allowindex = true)
|
||||
{
|
||||
double vlam[3];
|
||||
if(searchtree)
|
||||
const_cast<Mesh&>(mesh).BuildElementSearchTree(2);
|
||||
auto velement = Find2dElement(mesh, p, vlam, nullopt, searchtree ? mesh.GetSurfaceElementSearchTree() : nullptr, allowindex);
|
||||
if(!velement.IsValid())
|
||||
int velement = Find2dElement(mesh, p, vlam, NULL, searchtree, allowindex);
|
||||
if(velement == 0)
|
||||
return 0;
|
||||
|
||||
vlam[2] = 1.-vlam[0] - vlam[1];
|
||||
@ -233,8 +235,8 @@ namespace netgen
|
||||
: topology(*this), surfarea(*this)
|
||||
{
|
||||
lochfunc = {nullptr};
|
||||
for(auto i : Range(4))
|
||||
elementsearchtreets[i] = NextTimeStamp();
|
||||
elementsearchtree = nullptr;
|
||||
elementsearchtreets = NextTimeStamp();
|
||||
majortimestamp = timestamp = NextTimeStamp();
|
||||
hglob = 1e10;
|
||||
hmin = 0;
|
||||
@ -4819,21 +4821,6 @@ namespace netgen
|
||||
for (auto & seg : LineSegments())
|
||||
seg.si = seg.edgenr;
|
||||
}
|
||||
if (dimension == 3 && dim == 1)
|
||||
{
|
||||
for(auto str : materials)
|
||||
delete str;
|
||||
materials.SetSize(0);
|
||||
for(auto str : bcnames)
|
||||
delete str;
|
||||
bcnames.SetSize(0);
|
||||
for(auto str: cd2names)
|
||||
materials.Append(str);
|
||||
cd2names.SetSize(0);
|
||||
for(auto str : cd3names)
|
||||
bcnames.Append(str);
|
||||
cd3names.SetSize(0);
|
||||
}
|
||||
dimension = dim;
|
||||
}
|
||||
|
||||
@ -5264,95 +5251,122 @@ namespace netgen
|
||||
RebuildSurfaceElementLists();
|
||||
}
|
||||
|
||||
void Mesh :: BuildElementSearchTree (int dim)
|
||||
void Mesh :: BuildElementSearchTree ()
|
||||
{
|
||||
if(dim < 2)
|
||||
return;
|
||||
if (elementsearchtreets[dim] == GetTimeStamp())
|
||||
return;
|
||||
if (elementsearchtreets == GetTimeStamp()) return;
|
||||
|
||||
{
|
||||
std::lock_guard<std::mutex> guard(buildsearchtree_mutex);
|
||||
// check again to see if some other thread built while waiting for lock
|
||||
if (elementsearchtreets[dim] == GetTimeStamp()) return;
|
||||
|
||||
elementsearchtreets[dim] = GetTimeStamp();
|
||||
|
||||
PrintMessage (4, "Rebuild element searchtree dim " + ToString(dim));
|
||||
if (elementsearchtreets != GetTimeStamp())
|
||||
{
|
||||
NgLock lock(mutex);
|
||||
lock.Lock();
|
||||
|
||||
PrintMessage (4, "Rebuild element searchtree");
|
||||
|
||||
elementsearchtree = nullptr;
|
||||
|
||||
int ne = (dimension == 2) ? GetNSE() : GetNE();
|
||||
if (dimension == 3 && !GetNE() && GetNSE())
|
||||
ne = GetNSE();
|
||||
|
||||
Point3d pmin, pmax;
|
||||
GetBox(pmin, pmax);
|
||||
Box<3> box(pmin, pmax);
|
||||
if (dim == 3)
|
||||
elementsearchtree_vol = make_unique<BoxTree<3, ElementIndex>>(box);
|
||||
else
|
||||
elementsearchtree_surf = make_unique<BoxTree<3, SurfaceElementIndex>>(box);
|
||||
|
||||
if (dim == 3)
|
||||
{
|
||||
for(auto ei : volelements.Range())
|
||||
if (ne)
|
||||
{
|
||||
const auto& el = volelements[ei];
|
||||
Box<3> box (Box<3>::EMPTY_BOX);
|
||||
for (auto pi : el.PNums())
|
||||
box.Add (points[pi]);
|
||||
|
||||
if(el.IsCurved() && curvedelems->IsElementCurved(ei))
|
||||
if (dimension == 2 || (dimension == 3 && !GetNE()) )
|
||||
{
|
||||
// add edge/face midpoints to box
|
||||
auto eltype = el.GetType();
|
||||
const auto verts = topology.GetVertices(eltype);
|
||||
|
||||
const auto edges = FlatArray<const ELEMENT_EDGE>(topology.GetNEdges(eltype), topology.GetEdges0(eltype));
|
||||
for (const auto & edge: edges) {
|
||||
netgen::Point<3> lam = netgen::Point<3>(0.5* (verts[edge[0]] + verts[edge[1]]));
|
||||
auto p = netgen::Point<3>(0.0);
|
||||
curvedelems->CalcElementTransformation(lam,ei,p);
|
||||
box.Add(p);
|
||||
}
|
||||
|
||||
const auto faces = FlatArray<const ELEMENT_FACE>(topology.GetNFaces(eltype), topology.GetFaces0(eltype));
|
||||
for (const auto & face: faces) {
|
||||
netgen::Vec<3> lam = netgen::Vec<3>(verts[face[0]] + verts[face[1]] + verts[face[2]]);
|
||||
if(face[3] != -1) {
|
||||
lam += netgen::Vec<3>(verts[face[3]]);
|
||||
lam *= 0.25;
|
||||
}
|
||||
else
|
||||
lam *= 1.0/3;
|
||||
auto p = netgen::Point<3>(0.0);
|
||||
curvedelems->CalcElementTransformation(netgen::Point<3>(lam),ei,p);
|
||||
box.Add(p);
|
||||
}
|
||||
}
|
||||
box.Scale(1.2);
|
||||
elementsearchtree_vol -> Insert (box, ei);
|
||||
}
|
||||
}
|
||||
else if (dim == 2)
|
||||
{
|
||||
for (auto ei : Range(surfelements))
|
||||
{
|
||||
const auto& el = surfelements[ei];
|
||||
Box<3> box (Box<3>::EMPTY_BOX);
|
||||
for (auto pi : el.PNums())
|
||||
box.Add (points[pi]);
|
||||
|
||||
if(el.IsCurved() && curvedelems->IsSurfaceElementCurved(ei))
|
||||
{
|
||||
netgen::Point<2> lami [4] = {netgen::Point<2>(0.5,0), netgen::Point<2>(0,0.5), netgen::Point<2>(0.5,0.5), netgen::Point<2>(1./3,1./3)};
|
||||
for (auto lam : lami)
|
||||
Box<3> box (Box<3>::EMPTY_BOX);
|
||||
for (SurfaceElementIndex sei = 0; sei < ne; sei++)
|
||||
// box.Add (points[surfelements[sei].PNums()]);
|
||||
for (auto pi : surfelements[sei].PNums())
|
||||
box.Add (points[pi]);
|
||||
|
||||
box.Increase (1.01 * box.Diam());
|
||||
elementsearchtree = make_unique<BoxTree<3>> (box);
|
||||
|
||||
for (SurfaceElementIndex sei = 0; sei < ne; sei++)
|
||||
{
|
||||
netgen::Point<3> x;
|
||||
Mat<3,2> Jac;
|
||||
// box.Set (points[surfelements[sei].PNums()]);
|
||||
|
||||
curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
|
||||
box.Add (x);
|
||||
Box<3> box (Box<3>::EMPTY_BOX);
|
||||
for (auto pi : surfelements[sei].PNums())
|
||||
box.Add (points[pi]);
|
||||
|
||||
auto & el = surfelements[sei];
|
||||
if(el.IsCurved() && curvedelems->IsSurfaceElementCurved(sei))
|
||||
{
|
||||
netgen::Point<2> lami [4] = {netgen::Point<2>(0.5,0), netgen::Point<2>(0,0.5), netgen::Point<2>(0.5,0.5), netgen::Point<2>(1./3,1./3)};
|
||||
for (auto lam : lami)
|
||||
{
|
||||
netgen::Point<3> x;
|
||||
Mat<3,2> Jac;
|
||||
|
||||
curvedelems->CalcSurfaceTransformation(lam,sei,x,Jac);
|
||||
box.Add (x);
|
||||
}
|
||||
box.Scale(1.2);
|
||||
}
|
||||
elementsearchtree -> Insert (box, sei+1);
|
||||
}
|
||||
box.Scale(1.2);
|
||||
}
|
||||
elementsearchtree_surf -> Insert (box, ei);
|
||||
else
|
||||
{
|
||||
Box<3> box (Box<3>::EMPTY_BOX);
|
||||
for (ElementIndex ei = 0; ei < ne; ei++)
|
||||
// box.Add (points[volelements[ei].PNums()]);
|
||||
for (auto pi : volelements[ei].PNums())
|
||||
box.Add (points[pi]);
|
||||
|
||||
box.Increase (1.01 * box.Diam());
|
||||
elementsearchtree = make_unique<BoxTree<3>> (box);
|
||||
|
||||
for (ElementIndex ei = 0; ei < ne; ei++)
|
||||
{
|
||||
// box.Set (points[volelements[ei].PNums()]);
|
||||
|
||||
Box<3> box (Box<3>::EMPTY_BOX);
|
||||
for (auto pi : volelements[ei].PNums())
|
||||
box.Add (points[pi]);
|
||||
|
||||
auto & el = volelements[ei];
|
||||
if(el.IsCurved() && curvedelems->IsElementCurved(ei))
|
||||
{
|
||||
// add edge/face midpoints to box
|
||||
auto eltype = el.GetType();
|
||||
const auto verts = topology.GetVertices(eltype);
|
||||
|
||||
|
||||
const auto edges = FlatArray<const ELEMENT_EDGE>(topology.GetNEdges(eltype), topology.GetEdges0(eltype));
|
||||
for (const auto & edge: edges) {
|
||||
netgen::Point<3> lam = netgen::Point<3>(0.5* (verts[edge[0]] + verts[edge[1]]));
|
||||
auto p = netgen::Point<3>(0.0);
|
||||
curvedelems->CalcElementTransformation(lam,ei,p);
|
||||
box.Add(p);
|
||||
}
|
||||
|
||||
const auto faces = FlatArray<const ELEMENT_FACE>(topology.GetNFaces(eltype), topology.GetFaces0(eltype));
|
||||
for (const auto & face: faces) {
|
||||
netgen::Vec<3> lam = netgen::Vec<3>(verts[face[0]] + verts[face[1]] + verts[face[2]]);
|
||||
if(face[3] != -1) {
|
||||
lam += netgen::Vec<3>(verts[face[3]]);
|
||||
lam *= 0.25;
|
||||
}
|
||||
else
|
||||
lam *= 1.0/3;
|
||||
auto p = netgen::Point<3>(0.0);
|
||||
curvedelems->CalcElementTransformation(netgen::Point<3>(lam),ei,p);
|
||||
box.Add(p);
|
||||
}
|
||||
|
||||
|
||||
box.Scale(1.2);
|
||||
}
|
||||
|
||||
|
||||
elementsearchtree -> Insert (box, ei+1);
|
||||
}
|
||||
}
|
||||
|
||||
elementsearchtreets = GetTimeStamp();
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -5392,7 +5406,7 @@ namespace netgen
|
||||
|
||||
bool Mesh :: PointContainedIn2DElement(const Point3d & p,
|
||||
double lami[3],
|
||||
SurfaceElementIndex ei,
|
||||
const int element,
|
||||
bool consider3D) const
|
||||
{
|
||||
Vec3d col1, col2, col3;
|
||||
@ -5403,9 +5417,9 @@ namespace netgen
|
||||
|
||||
|
||||
//SZ
|
||||
if(surfelements[ei].GetType()==QUAD)
|
||||
if(SurfaceElement(element).GetType()==QUAD)
|
||||
{
|
||||
const Element2d & el = surfelements[ei];
|
||||
const Element2d & el = SurfaceElement(element);
|
||||
|
||||
const Point3d & p1 = Point(el.PNum(1));
|
||||
const Point3d & p2 = Point(el.PNum(2));
|
||||
@ -5424,7 +5438,7 @@ namespace netgen
|
||||
int i = 0;
|
||||
while(delta > 1e-16 && i < maxits)
|
||||
{
|
||||
curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
|
||||
curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac);
|
||||
rhs = p - x;
|
||||
Jac.Solve(rhs,deltalam);
|
||||
lam += deltalam;
|
||||
@ -5753,7 +5767,7 @@ namespace netgen
|
||||
{
|
||||
// SurfaceElement(element).GetTets (loctets);
|
||||
loctrigs.SetSize(1);
|
||||
loctrigs.Elem(1) = surfelements[ei];
|
||||
loctrigs.Elem(1) = SurfaceElement(element);
|
||||
|
||||
|
||||
|
||||
@ -5788,11 +5802,11 @@ namespace netgen
|
||||
//(*testout) << "col1 " << col1 << " col2 " << col2 << " col3 " << col3 << " rhs " << rhs << endl;
|
||||
//(*testout) << "sol " << sol << endl;
|
||||
|
||||
if (surfelements[ei].GetType() ==TRIG6 || curvedelems->IsSurfaceElementCurved(ei))
|
||||
if (SurfaceElement(element).GetType() ==TRIG6 || curvedelems->IsSurfaceElementCurved(element-1))
|
||||
{
|
||||
// netgen::Point<2> lam(1./3,1./3);
|
||||
netgen::Point<2> lam(sol.X(), sol.Y());
|
||||
if(surfelements[ei].GetType() != TRIG6)
|
||||
if(SurfaceElement(element).GetType() != TRIG6)
|
||||
{
|
||||
lam[0] = 1-sol.X()-sol.Y();
|
||||
lam[1] = sol.X();
|
||||
@ -5811,7 +5825,7 @@ namespace netgen
|
||||
const int maxits = 30;
|
||||
while(delta > 1e-16 && i<maxits)
|
||||
{
|
||||
curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
|
||||
curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac);
|
||||
rhs = p-x;
|
||||
Jac.Solve(rhs,deltalam);
|
||||
|
||||
@ -5830,7 +5844,7 @@ namespace netgen
|
||||
sol.X() = lam(0);
|
||||
sol.Y() = lam(1);
|
||||
|
||||
if (surfelements[ei].GetType() !=TRIG6 )
|
||||
if (SurfaceElement(element).GetType() !=TRIG6 )
|
||||
{
|
||||
sol.Z() = sol.X();
|
||||
sol.X() = sol.Y();
|
||||
@ -5862,7 +5876,7 @@ namespace netgen
|
||||
|
||||
bool Mesh :: PointContainedIn3DElement(const Point3d & p,
|
||||
double lami[3],
|
||||
ElementIndex ei) const
|
||||
const int element) const
|
||||
{
|
||||
//bool oldresult = PointContainedIn3DElementOld(p,lami,element);
|
||||
//(*testout) << "old result: " << oldresult
|
||||
@ -5873,7 +5887,7 @@ namespace netgen
|
||||
|
||||
|
||||
const double eps = 1.e-4;
|
||||
const Element & el = volelements[ei];
|
||||
const Element & el = VolumeElement(element);
|
||||
|
||||
netgen::Point<3> lam = 0.0;
|
||||
|
||||
@ -5908,7 +5922,7 @@ namespace netgen
|
||||
const int maxits = 30;
|
||||
while(delta > 1e-16 && i<maxits)
|
||||
{
|
||||
curvedelems->CalcElementTransformation(lam,ei,x,Jac);
|
||||
curvedelems->CalcElementTransformation(lam,element-1,x,Jac);
|
||||
rhs = p-x;
|
||||
Jac.Solve(rhs,deltalam);
|
||||
|
||||
@ -6030,72 +6044,86 @@ namespace netgen
|
||||
}
|
||||
|
||||
|
||||
ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
|
||||
double* lami,
|
||||
bool build_searchtree,
|
||||
int index,
|
||||
bool allowindex) const
|
||||
int Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
|
||||
double lami[3],
|
||||
bool build_searchtree,
|
||||
const int index,
|
||||
const bool allowindex) const
|
||||
{
|
||||
if(index != -1)
|
||||
{
|
||||
Array<int> dummy(1);
|
||||
NgArray<int> dummy(1);
|
||||
dummy[0] = index;
|
||||
return GetElementOfPoint(p,lami,dummy,build_searchtree,allowindex);
|
||||
return GetElementOfPoint(p,lami,&dummy,build_searchtree,allowindex);
|
||||
}
|
||||
else
|
||||
return GetElementOfPoint(p,lami,nullopt,build_searchtree,allowindex);
|
||||
return GetElementOfPoint(p,lami,NULL,build_searchtree,allowindex);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
|
||||
double* lami,
|
||||
std::optional<FlatArray<int>> indices,
|
||||
bool build_searchtree,
|
||||
bool allowindex) const
|
||||
int Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
|
||||
double lami[3],
|
||||
const NgArray<int> * const indices,
|
||||
bool build_searchtree,
|
||||
const bool allowindex) const
|
||||
{
|
||||
if ( (dimension == 2 && !GetNSE()) ||
|
||||
(dimension == 3 && !GetNE() && !GetNSE()) )
|
||||
return -1;
|
||||
|
||||
if (build_searchtree)
|
||||
const_cast<Mesh&>(*this).BuildElementSearchTree (3);
|
||||
return Find3dElement(*this, p, lami, indices, elementsearchtree_vol.get(), allowindex);
|
||||
const_cast<Mesh&>(*this).BuildElementSearchTree ();
|
||||
|
||||
if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE()))
|
||||
return Find2dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
|
||||
|
||||
return Find3dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
|
||||
}
|
||||
|
||||
|
||||
|
||||
SurfaceElementIndex Mesh ::
|
||||
GetSurfaceElementOfPoint (const netgen::Point<3> & p,
|
||||
double* lami,
|
||||
bool build_searchtree,
|
||||
int index,
|
||||
bool allowindex) const
|
||||
int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p,
|
||||
double lami[3],
|
||||
bool build_searchtree,
|
||||
const int index,
|
||||
const bool allowindex) const
|
||||
{
|
||||
if(index != -1)
|
||||
if(index != -1)
|
||||
{
|
||||
Array<int> dummy(1);
|
||||
NgArray<int> dummy(1);
|
||||
dummy[0] = index;
|
||||
return GetSurfaceElementOfPoint(p,lami,dummy,build_searchtree,allowindex);
|
||||
return GetSurfaceElementOfPoint(p,lami,&dummy,build_searchtree,allowindex);
|
||||
}
|
||||
else
|
||||
return GetSurfaceElementOfPoint(p,lami,nullopt,build_searchtree,allowindex);
|
||||
return GetSurfaceElementOfPoint(p,lami,NULL,build_searchtree,allowindex);
|
||||
}
|
||||
|
||||
SurfaceElementIndex Mesh ::
|
||||
GetSurfaceElementOfPoint (const netgen::Point<3> & p,
|
||||
double* lami,
|
||||
std::optional<FlatArray<int>> indices,
|
||||
bool build_searchtree,
|
||||
bool allowindex) const
|
||||
|
||||
|
||||
|
||||
int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p,
|
||||
double lami[3],
|
||||
const NgArray<int> * const indices,
|
||||
bool build_searchtree,
|
||||
const bool allowindex) const
|
||||
{
|
||||
if (build_searchtree)
|
||||
const_cast<Mesh&>(*this).BuildElementSearchTree(2);
|
||||
return Find2dElement(*this, p, lami, indices, elementsearchtree_surf.get(), allowindex);
|
||||
if (!GetNE() && build_searchtree)
|
||||
const_cast<Mesh&>(*this).BuildElementSearchTree ();
|
||||
|
||||
if (dimension == 2)
|
||||
return Find1dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
|
||||
else
|
||||
return Find2dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2,
|
||||
Array<ElementIndex> & locels) const
|
||||
NgArray<int> & locels) const
|
||||
{
|
||||
elementsearchtree_vol->GetIntersecting (p1, p2, locels);
|
||||
elementsearchtree->GetIntersecting (p1, p2, locels);
|
||||
}
|
||||
|
||||
void Mesh :: SplitIntoParts()
|
||||
@ -7413,10 +7441,9 @@ namespace netgen
|
||||
}
|
||||
|
||||
string Mesh :: defaultmat = "default";
|
||||
string_view Mesh :: defaultmat_sv = "default";
|
||||
const string & Mesh :: GetMaterial (int domnr) const
|
||||
{
|
||||
if (domnr <= materials.Size() && materials[domnr-1])
|
||||
if (domnr <= materials.Size())
|
||||
return *materials[domnr-1];
|
||||
static string emptystring("default");
|
||||
return emptystring;
|
||||
|
||||
@ -128,13 +128,6 @@ namespace netgen
|
||||
*/
|
||||
NgArray<EdgeDescriptor> edgedecoding;
|
||||
|
||||
Array<string*> region_name_cd[4];
|
||||
Array<string*> & materials = region_name_cd[0];
|
||||
Array<string*> & bcnames = region_name_cd[1];
|
||||
Array<string*> & cd2names = region_name_cd[2];
|
||||
Array<string*> & cd3names = region_name_cd[3];
|
||||
|
||||
/*
|
||||
/// sub-domain materials
|
||||
Array<string*> materials;
|
||||
|
||||
@ -146,8 +139,7 @@ namespace netgen
|
||||
|
||||
/// labels for co dim 3 bbboundary conditions
|
||||
Array<string*> cd3names;
|
||||
*/
|
||||
|
||||
|
||||
/// Periodic surface, close surface, etc. identifications
|
||||
unique_ptr<Identifications> ident;
|
||||
|
||||
@ -156,10 +148,9 @@ namespace netgen
|
||||
int numvertices;
|
||||
|
||||
/// geometric search tree for interval intersection search
|
||||
unique_ptr<BoxTree<3, ElementIndex>> elementsearchtree_vol;
|
||||
unique_ptr<BoxTree<3, SurfaceElementIndex>> elementsearchtree_surf;
|
||||
unique_ptr<BoxTree<3>> elementsearchtree;
|
||||
/// time stamp for tree
|
||||
mutable size_t elementsearchtreets[4];
|
||||
mutable int elementsearchtreets;
|
||||
|
||||
/// element -> face, element -> edge etc ...
|
||||
MeshTopology topology;
|
||||
@ -209,11 +200,11 @@ namespace netgen
|
||||
|
||||
DLL_HEADER bool PointContainedIn2DElement(const Point3d & p,
|
||||
double lami[3],
|
||||
SurfaceElementIndex element,
|
||||
const int element,
|
||||
bool consider3D = false) const;
|
||||
DLL_HEADER bool PointContainedIn3DElement(const Point3d & p,
|
||||
double lami[3],
|
||||
ElementIndex element) const;
|
||||
const int element) const;
|
||||
DLL_HEADER bool PointContainedIn3DElementOld(const Point3d & p,
|
||||
double lami[3],
|
||||
const int element) const;
|
||||
@ -672,48 +663,35 @@ namespace netgen
|
||||
|
||||
|
||||
/// build box-search tree
|
||||
DLL_HEADER void BuildElementSearchTree (int dim);
|
||||
BoxTree<3, ElementIndex>* GetElementSearchTree () const
|
||||
{
|
||||
return elementsearchtree_vol.get();
|
||||
}
|
||||
|
||||
BoxTree<3, SurfaceElementIndex>* GetSurfaceElementSearchTree () const
|
||||
{
|
||||
return elementsearchtree_surf.get();
|
||||
}
|
||||
DLL_HEADER void BuildElementSearchTree ();
|
||||
|
||||
void SetPointSearchStartElement(const int el) const {ps_startelement = el;}
|
||||
|
||||
/// gives element of point, barycentric coordinates
|
||||
DLL_HEADER ElementIndex
|
||||
GetElementOfPoint (const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
bool build_searchtree = false,
|
||||
int index = -1,
|
||||
bool allowindex = true) const;
|
||||
DLL_HEADER ElementIndex
|
||||
GetElementOfPoint (const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
std::optional<FlatArray<int>> indices,
|
||||
bool build_searchtree = 0,
|
||||
bool allowindex = true) const;
|
||||
DLL_HEADER SurfaceElementIndex
|
||||
GetSurfaceElementOfPoint (const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
bool build_searchtree = false,
|
||||
int index = -1,
|
||||
bool allowindex = true) const;
|
||||
DLL_HEADER SurfaceElementIndex
|
||||
GetSurfaceElementOfPoint (const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
std::optional<FlatArray<int>> indices,
|
||||
bool build_searchtree = false,
|
||||
bool allowindex = true) const;
|
||||
DLL_HEADER int GetElementOfPoint (const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
bool build_searchtree = 0,
|
||||
const int index = -1,
|
||||
const bool allowindex = true) const;
|
||||
DLL_HEADER int GetElementOfPoint (const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
const NgArray<int> * const indices,
|
||||
bool build_searchtree = 0,
|
||||
const bool allowindex = true) const;
|
||||
DLL_HEADER int GetSurfaceElementOfPoint (const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
bool build_searchtree = 0,
|
||||
const int index = -1,
|
||||
const bool allowindex = true) const;
|
||||
DLL_HEADER int GetSurfaceElementOfPoint (const netgen::Point<3> & p,
|
||||
double * lami,
|
||||
const NgArray<int> * const indices,
|
||||
bool build_searchtree = 0,
|
||||
const bool allowindex = true) const;
|
||||
|
||||
/// give list of vol elements which are int the box(p1,p2)
|
||||
void GetIntersectingVolEls(const Point3d& p1, const Point3d& p2,
|
||||
Array<ElementIndex> & locels) const;
|
||||
NgArray<int> & locels) const;
|
||||
|
||||
///
|
||||
int AddFaceDescriptor(const FaceDescriptor& fd)
|
||||
@ -783,15 +761,6 @@ namespace netgen
|
||||
std::string_view GetRegionName(SegmentIndex ei) const { return GetRegionName((*this)[ei]); }
|
||||
std::string_view GetRegionName(SurfaceElementIndex ei) const { return GetRegionName((*this)[ei]); }
|
||||
std::string_view GetRegionName(ElementIndex ei) const { return GetRegionName((*this)[ei]); }
|
||||
|
||||
DLL_HEADER static string_view defaultmat_sv;
|
||||
std::string_view GetRegionName (int dim, int domnr) // 1-based domnr
|
||||
{
|
||||
domnr--;
|
||||
auto & names = region_name_cd[dimension-dim];
|
||||
if (domnr < names.Size() && names[domnr]) return *names[domnr];
|
||||
return defaultmat_sv;
|
||||
}
|
||||
|
||||
///
|
||||
void ClearFaceDescriptors()
|
||||
@ -1067,6 +1036,7 @@ namespace netgen
|
||||
}
|
||||
|
||||
DLL_HEADER void AddFacesBetweenDomains(Mesh & mesh);
|
||||
|
||||
}
|
||||
|
||||
#endif // NETGEN_MESHCLASS_HPP
|
||||
|
||||
@ -783,7 +783,7 @@ namespace netgen
|
||||
free_segs.Append(segi);
|
||||
|
||||
auto get_nonconforming = [&] (const auto & p2el) {
|
||||
Array<SegmentIndex> nonconforming;
|
||||
Array<size_t> nonconforming;
|
||||
|
||||
for (auto segi : free_segs) {
|
||||
auto seg = mesh[segi];
|
||||
@ -805,56 +805,6 @@ namespace netgen
|
||||
auto split_segment = [&] (SegmentIndex segi, const auto & p2el) {
|
||||
auto seg = mesh[segi];
|
||||
auto p_new = Center(mesh[seg[0]], mesh[seg[1]]);
|
||||
double lam[3];
|
||||
ElementIndex ei_start = mesh.GetElementOfPoint(p_new, lam, false, domain);
|
||||
|
||||
if(!ei_start.IsValid()) {
|
||||
PrintMessage(1, "Could not find volume element with new point");
|
||||
return;
|
||||
}
|
||||
|
||||
if(mesh[ei_start].IsDeleted())
|
||||
return;
|
||||
|
||||
double max_inside = -1.;
|
||||
ElementIndex ei_max_inside = ElementIndex::INVALID;
|
||||
|
||||
// search for adjacent volume element, where the new point is "most inside",
|
||||
// i.e. the minimal barycentric coordinate is maximal
|
||||
for(auto pi : mesh[ei_start].PNums()) {
|
||||
for(auto ei1 : p2el[pi]) {
|
||||
double lam[3];
|
||||
|
||||
if(mesh[ei1].IsDeleted())
|
||||
return;
|
||||
if(!mesh.PointContainedIn3DElement(p_new, lam, ei1))
|
||||
continue;
|
||||
|
||||
double inside = min(min(lam[0], lam[1]), min(lam[2], 1.0-lam[0]-lam[1]));
|
||||
if(inside > max_inside) {
|
||||
max_inside = inside;
|
||||
ei_max_inside = ei1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if(max_inside < 1e-4) {
|
||||
PrintMessage(3, "Could not find volume element with new point inside");
|
||||
return;
|
||||
}
|
||||
|
||||
// split tet into 4 new tests, with new point inside
|
||||
auto el = mesh[ei_max_inside];
|
||||
if(el.GetNP() != 4) {
|
||||
PrintMessage(3, "Only tet elements are supported to split around free segments");
|
||||
return;
|
||||
}
|
||||
|
||||
if(el.IsDeleted()) {
|
||||
PrintMessage(3,"Element to split is already deleted");
|
||||
return;
|
||||
}
|
||||
|
||||
auto pi_new = mesh.AddPoint(p_new);
|
||||
auto seg_new0 = seg;
|
||||
auto seg_new1 = seg;
|
||||
@ -865,6 +815,45 @@ namespace netgen
|
||||
mesh.AddSegment(seg_new0);
|
||||
mesh.AddSegment(seg_new1);
|
||||
|
||||
double lam[3];
|
||||
ElementIndex ei_start = mesh.GetElementOfPoint(p_new, lam, false, domain);
|
||||
|
||||
if(ei_start == 0) {
|
||||
PrintMessage(1, "Could not find volume element with new point");
|
||||
return;
|
||||
}
|
||||
ei_start -= 1;
|
||||
|
||||
double max_inside = -1.;
|
||||
ElementIndex ei_max_inside = -1;
|
||||
|
||||
// search for adjacent volume element, where the new point is "most inside",
|
||||
// i.e. the minimal barycentric coordinate is maximal
|
||||
for(auto pi : mesh[ei_start].PNums()) {
|
||||
for(auto ei1 : p2el[pi]) {
|
||||
double lam[3];
|
||||
if(!mesh.PointContainedIn3DElement(p_new, lam, ei1+1))
|
||||
continue;
|
||||
|
||||
double inside = min(min(lam[0], lam[1]), min(lam[2], 1.0-lam[0]-lam[1]));
|
||||
if(inside > max_inside) {
|
||||
max_inside = inside;
|
||||
ei_max_inside = ei1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// split tet into 4 new tests, with new point inside
|
||||
auto el = mesh[ei_max_inside];
|
||||
if(el.GetNP() != 4) {
|
||||
PrintMessage(1, "Only tet elements are supported to split around free segments");
|
||||
return;
|
||||
}
|
||||
|
||||
if(el.IsDeleted()) {
|
||||
PrintMessage(1,"Element to split is already deleted");
|
||||
return;
|
||||
}
|
||||
|
||||
int pmap[4][4] = {
|
||||
{0,1,2,4},
|
||||
@ -908,7 +897,7 @@ namespace netgen
|
||||
|
||||
for ([[maybe_unused]] auto i : Range(3)) {
|
||||
optmesh.ImproveMesh();
|
||||
optmesh.SwapImprove2(true);
|
||||
optmesh.SwapImprove2 ();
|
||||
optmesh.ImproveMesh();
|
||||
optmesh.SwapImprove();
|
||||
optmesh.ImproveMesh();
|
||||
@ -921,7 +910,7 @@ namespace netgen
|
||||
auto bad_segs = get_nonconforming(p2el);
|
||||
|
||||
if(bad_segs.Size() > 0) {
|
||||
auto bad_seg = mesh[bad_segs[0]];
|
||||
auto bad_seg = mesh[free_segs[bad_segs[0]]];
|
||||
if(debugparam.write_mesh_on_error)
|
||||
mesh.Save("free_segment_not_conformed_dom_"+ToString(domain)+"_seg_"+ToString(bad_seg[0])+"_"+ToString(bad_seg[1])+".vol.gz");
|
||||
throw Exception("Segment not resolved in volume mesh in domain " + ToString(domain)+ ", seg: " + ToString(bad_seg));
|
||||
|
||||
@ -1635,7 +1635,7 @@ namespace netgen
|
||||
Point<3> pnt;
|
||||
double h;
|
||||
int layer = 1;
|
||||
MeshSizePoint (Point<3> pnt_, double h_, int layer_ = 1) : pnt(pnt_), h(h_), layer(layer_) { ; }
|
||||
MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; }
|
||||
MeshSizePoint () = default;
|
||||
MeshSizePoint (const MeshSizePoint &) = default;
|
||||
MeshSizePoint (MeshSizePoint &&) = default;
|
||||
|
||||
@ -966,7 +966,7 @@ namespace netgen
|
||||
self.ident = make_unique<Identifications> (self);
|
||||
self.topology = MeshTopology(*this);
|
||||
self.topology.Update();
|
||||
self.BuildElementSearchTree(3);
|
||||
self.BuildElementSearchTree();
|
||||
|
||||
// const_cast<Mesh&>(*this).DeleteMesh();
|
||||
|
||||
|
||||
@ -1300,10 +1300,6 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
|
||||
.def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries,
|
||||
py::arg("identification_name"), py::arg("face1"), py::arg("mapping"),
|
||||
py::arg("point_tolerance") = -1.)
|
||||
.def("GetCurveOrder", [] (Mesh & self)
|
||||
{
|
||||
return self.GetCurvedElements().GetOrder();
|
||||
})
|
||||
.def("GetNrIdentifications", [](Mesh& self)
|
||||
{
|
||||
return self.GetIdentifications().GetMaxNr();
|
||||
@ -1484,8 +1480,7 @@ py::arg("point_tolerance") = -1.)
|
||||
}))
|
||||
*/
|
||||
|
||||
.def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>(),
|
||||
py::arg("dim")=3)
|
||||
.def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>())
|
||||
|
||||
.def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array<int>{})
|
||||
.def ("BoundaryLayer", [](Mesh & self, variant<string, int, std::vector<int>> boundary,
|
||||
@ -1687,25 +1682,25 @@ py::arg("point_tolerance") = -1.)
|
||||
return mp;
|
||||
}), py::arg("mp")=nullptr, meshingparameter_description.c_str())
|
||||
.def("__str__", &ToString<MP>)
|
||||
.def("RestrictH", [](MP & mp, double x, double y, double z, double h, int layer)
|
||||
.def("RestrictH", [](MP & mp, double x, double y, double z, double h)
|
||||
{
|
||||
mp.meshsize_points.Append ( MeshingParameters::MeshSizePoint(Point<3> (x,y,z), h, layer));
|
||||
}, py::arg("x"), py::arg("y"), py::arg("z"), py::arg("h"), py::arg("layer")=1
|
||||
mp.meshsize_points.Append ( MeshingParameters::MeshSizePoint(Point<3> (x,y,z), h));
|
||||
}, py::arg("x"), py::arg("y"), py::arg("z"), py::arg("h")
|
||||
)
|
||||
.def("RestrictH", [](MP & mp, const Point<3>& p, double h, int layer)
|
||||
.def("RestrictH", [](MP & mp, const Point<3>& p, double h)
|
||||
{
|
||||
mp.meshsize_points.Append ({p, h, layer});
|
||||
}, py::arg("p"), py::arg("h"), py::arg("layer")=1)
|
||||
mp.meshsize_points.Append ({p, h});
|
||||
}, py::arg("p"), py::arg("h"))
|
||||
.def("RestrictHLine", [](MP& mp, const Point<3>& p1, const Point<3>& p2,
|
||||
double maxh, int layer)
|
||||
double maxh)
|
||||
{
|
||||
int steps = int(Dist(p1, p2) / maxh) + 2;
|
||||
auto v = p2 - p1;
|
||||
for (int i = 0; i <= steps; i++)
|
||||
{
|
||||
mp.meshsize_points.Append({p1 + double(i)/steps * v, maxh, layer});
|
||||
mp.meshsize_points.Append({p1 + double(i)/steps * v, maxh});
|
||||
}
|
||||
}, py::arg("p1"), py::arg("p2"), py::arg("maxh"), py::arg("layer")=1)
|
||||
}, py::arg("p1"), py::arg("p2"), py::arg("maxh"))
|
||||
;
|
||||
|
||||
m.def("SetTestoutFile", FunctionPointer ([] (const string & filename)
|
||||
|
||||
@ -2288,22 +2288,10 @@ namespace netgen
|
||||
XCAFPrs::CollectStyleSettings(label, loc, set);
|
||||
XCAFPrs_Style aStyle;
|
||||
set.FindFromKey(e.Current(), aStyle);
|
||||
|
||||
auto & prop = OCCGeometry::GetProperties(e.Current());
|
||||
if(aStyle.IsSetColorSurf())
|
||||
{
|
||||
for(TopExp_Explorer e2(e.Current(), TopAbs_FACE); e2.More(); e2.Next())
|
||||
{
|
||||
auto & prop = OCCGeometry::GetProperties(e2.Current());
|
||||
prop.col = step_utils::ReadColor(aStyle.GetColorSurfRGBA());
|
||||
}
|
||||
}
|
||||
if(aStyle.IsSetColorCurv())
|
||||
{
|
||||
for(TopExp_Explorer e2(e.Current(), TopAbs_EDGE); e2.More(); e2.Next())
|
||||
{
|
||||
auto & prop = OCCGeometry::GetProperties(e2.Current());
|
||||
prop.col = step_utils::ReadColor(aStyle.GetColorSurfRGBA());
|
||||
}
|
||||
}
|
||||
prop.col = step_utils::ReadColor(aStyle.GetColorSurfRGBA());
|
||||
}
|
||||
|
||||
// load names
|
||||
|
||||
@ -21,9 +21,6 @@
|
||||
#include <XCAFDoc_DocumentTool.hxx>
|
||||
#include <XCAFDoc_MaterialTool.hxx>
|
||||
#include <XCAFDoc_ShapeTool.hxx>
|
||||
#include <TopoDS_Edge.hxx>
|
||||
#include <BRepAdaptor_Curve.hxx>
|
||||
#include <GCPnts_TangentialDeflection.hxx>
|
||||
|
||||
using namespace netgen;
|
||||
|
||||
@ -168,40 +165,10 @@ DLL_HEADER void ExportNgOCC(py::module &m)
|
||||
{
|
||||
ng_geometry = geo;
|
||||
})
|
||||
.def_property_readonly("solids", [](shared_ptr<OCCGeometry> geo)
|
||||
{
|
||||
ListOfShapes solids;
|
||||
for (int i = 1; i <= geo->somap.Extent(); i++)
|
||||
solids.push_back(geo->somap(i));
|
||||
return solids;
|
||||
}, "Get solids in order that they will be in the mesh")
|
||||
.def_property_readonly("faces", [](shared_ptr<OCCGeometry> geo)
|
||||
{
|
||||
ListOfShapes faces;
|
||||
for (int i = 1; i <= geo->fmap.Extent(); i++)
|
||||
faces.push_back(geo->fmap(i));
|
||||
return faces;
|
||||
}, "Get faces in order that they will be in the mesh")
|
||||
.def_property_readonly("edges", [](shared_ptr<OCCGeometry> geo)
|
||||
{
|
||||
ListOfShapes edges;
|
||||
for (int i = 1; i <= geo->emap.Extent(); i++)
|
||||
edges.push_back(geo->emap(i));
|
||||
return edges;
|
||||
}, "Get edges in order that they will be in the mesh")
|
||||
.def_property_readonly("vertices", [](shared_ptr<OCCGeometry> geo)
|
||||
{
|
||||
ListOfShapes vertices;
|
||||
for (int i = 1; i <= geo->vmap.Extent(); i++)
|
||||
vertices.push_back(geo->vmap(i));
|
||||
return vertices;
|
||||
}, "Get vertices in order that they will be in the mesh")
|
||||
.def("_visualizationData", [] (shared_ptr<OCCGeometry> occ_geo)
|
||||
{
|
||||
std::vector<float> vertices;
|
||||
std::vector<uint32_t> indices;
|
||||
std::vector<float> edges;
|
||||
std::vector<uint32_t> edge_indices;
|
||||
std::vector<int> trigs;
|
||||
std::vector<float> normals;
|
||||
std::vector<float> min = {std::numeric_limits<float>::max(),
|
||||
std::numeric_limits<float>::max(),
|
||||
@ -209,8 +176,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
|
||||
std::vector<float> max = {std::numeric_limits<float>::lowest(),
|
||||
std::numeric_limits<float>::lowest(),
|
||||
std::numeric_limits<float>::lowest()};
|
||||
std::vector<float> face_colors;
|
||||
std::vector<float> edge_colors;
|
||||
std::vector<string> surfnames;
|
||||
auto box = occ_geo->GetBoundingBox();
|
||||
for(int i = 0; i < 3; i++)
|
||||
{
|
||||
@ -222,67 +188,11 @@ DLL_HEADER void ExportNgOCC(py::module &m)
|
||||
gp_Pnt pnt;
|
||||
gp_Vec n;
|
||||
gp_Pnt p[3];
|
||||
for(int edge_index = 1; edge_index <= occ_geo->emap.Extent();
|
||||
edge_index++)
|
||||
{
|
||||
auto edge = TopoDS::Edge(occ_geo->emap(edge_index));
|
||||
if(OCCGeometry::HaveProperties(edge))
|
||||
{
|
||||
const auto& props = OCCGeometry::GetProperties(edge);
|
||||
if(props.col)
|
||||
edge_colors.insert(edge_colors.end(),
|
||||
{float((*props.col)[0]),
|
||||
float((*props.col)[1]),
|
||||
float((*props.col)[2]),
|
||||
float((*props.col)[3])});
|
||||
else
|
||||
edge_colors.insert(edge_colors.end(),{0.f,0.f,0.f,1.f});
|
||||
}
|
||||
else
|
||||
{
|
||||
edge_colors.insert(edge_colors.end(),{0.f,0.f,0.f,1.f});
|
||||
}
|
||||
BRepAdaptor_Curve adapt_crv = BRepAdaptor_Curve(edge);
|
||||
GCPnts_TangentialDeflection discretizer;
|
||||
discretizer.Initialize(adapt_crv, 0.09, 0.01);
|
||||
if (discretizer.NbPoints() > 1)
|
||||
{
|
||||
for (int j = 1; j <= discretizer.NbPoints()-1; ++j)
|
||||
{
|
||||
gp_Pnt p_0 = discretizer.Value(j);
|
||||
gp_Pnt p_1 = discretizer.Value(j+1);
|
||||
edges.insert(edges.end(),
|
||||
{float(p_0.X()),
|
||||
float(p_0.Y()),
|
||||
float(p_0.Z()),
|
||||
float(p_1.X()),
|
||||
float(p_1.Y()),
|
||||
float(p_1.Z())});
|
||||
edge_indices.push_back(uint32_t(edge_index-1));
|
||||
}
|
||||
}
|
||||
}
|
||||
int count = 0;
|
||||
for (int i = 1; i <= occ_geo->fmap.Extent(); i++)
|
||||
{
|
||||
surfnames.push_back("occ_surface" + to_string(i));
|
||||
auto face = TopoDS::Face(occ_geo->fmap(i));
|
||||
if (OCCGeometry::HaveProperties(face))
|
||||
{
|
||||
const auto& props = OCCGeometry::GetProperties(face);
|
||||
if(props.col)
|
||||
face_colors.insert(face_colors.end(),
|
||||
{float((*props.col)[0]),
|
||||
float((*props.col)[1]),
|
||||
float((*props.col)[2]),
|
||||
float((*props.col)[3])});
|
||||
else
|
||||
{
|
||||
face_colors.insert(face_colors.end(),{0.7,0.7,0.7,1.});
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
face_colors.insert(face_colors.end(),{0.7,0.7,0.7,1.});
|
||||
}
|
||||
auto surf = BRep_Tool::Surface(face);
|
||||
TopLoc_Location loc;
|
||||
BRepAdaptor_Surface sf(face, Standard_False);
|
||||
@ -290,7 +200,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
|
||||
Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
|
||||
if (triangulation.IsNull())
|
||||
cout << "cannot visualize face " << i << endl;
|
||||
indices.reserve(indices.size() + triangulation->NbTriangles());
|
||||
trigs.reserve(trigs.size() + triangulation->NbTriangles()*4);
|
||||
vertices.reserve(vertices.size() + triangulation->NbTriangles()*3*3);
|
||||
normals.reserve(normals.size() + triangulation->NbTriangles()*3*3);
|
||||
for (int j = 1; j < triangulation->NbTriangles()+1; j++)
|
||||
@ -298,13 +208,11 @@ DLL_HEADER void ExportNgOCC(py::module &m)
|
||||
auto triangle = triangulation->Triangle(j);
|
||||
for (int k = 1; k < 4; k++)
|
||||
p[k-1] = triangulation->Node(triangle(k)).Transformed(loc);
|
||||
indices.push_back(uint32_t(i-1));
|
||||
for (int k = 1; k < 4; k++)
|
||||
{
|
||||
vertices.insert(vertices.end(),{
|
||||
float(p[k-1].X()),
|
||||
float(p[k-1].Y()),
|
||||
float(p[k-1].Z())});
|
||||
vertices.insert(vertices.end(),{float(p[k-1].X()), float(p[k-1].Y()), float(p[k-1].Z())});
|
||||
trigs.insert(trigs.end(),{count, count+1, count+2,i});
|
||||
count += 3;
|
||||
uv = triangulation->UVNode(triangle(k));
|
||||
prop.SetParameters(uv.X(), uv.Y());
|
||||
if (prop.IsNormalDefined())
|
||||
@ -323,13 +231,12 @@ DLL_HEADER void ExportNgOCC(py::module &m)
|
||||
py::gil_scoped_acquire ac;
|
||||
py::dict res;
|
||||
py::list snames;
|
||||
for(auto name : surfnames)
|
||||
snames.append(py::cast(name));
|
||||
res["vertices"] = MoveToNumpy(vertices);
|
||||
res["edges"] = MoveToNumpy(edges);
|
||||
res["edge_indices"] = MoveToNumpy(edge_indices);
|
||||
res["edge_colors"] = MoveToNumpy(edge_colors);
|
||||
res["indices"] = MoveToNumpy(indices);
|
||||
res["triangles"] = MoveToNumpy(trigs);
|
||||
res["normals"] = MoveToNumpy(normals);
|
||||
res["face_colors"] = MoveToNumpy(face_colors);
|
||||
res["surfnames"] = snames;
|
||||
res["min"] = MoveToNumpy(min);
|
||||
res["max"] = MoveToNumpy(max);
|
||||
return res;
|
||||
|
||||
@ -14,9 +14,6 @@
|
||||
#pragma clang diagnostic push
|
||||
#pragma clang diagnostic ignored "-Wdeprecated-declarations"
|
||||
|
||||
#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
|
||||
#include <BinTools_ShapeWriter.hxx>
|
||||
#endif // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
|
||||
#include <BOPAlgo_Builder.hxx>
|
||||
#include <BOPTools_AlgoTools.hxx>
|
||||
#include <BRepAlgoAPI_Common.hxx>
|
||||
@ -52,7 +49,6 @@
|
||||
#include <BRepTools.hxx>
|
||||
#include <GCE2d_MakeArcOfCircle.hxx>
|
||||
#include <GCE2d_MakeCircle.hxx>
|
||||
#include <GCE2d_MakeEllipse.hxx>
|
||||
#include <GCE2d_MakeSegment.hxx>
|
||||
#include <GC_MakeArcOfCircle.hxx>
|
||||
#include <GC_MakeCircle.hxx>
|
||||
@ -602,19 +598,6 @@ public:
|
||||
return shared_from_this();
|
||||
}
|
||||
|
||||
auto Ellipse(double major, double minor)
|
||||
{
|
||||
Handle(Geom2d_Ellipse) ell_curve = GCE2d_MakeEllipse(localpos, major, minor).Value();
|
||||
|
||||
auto edge = BRepBuilderAPI_MakeEdge(ell_curve, surf).Edge();
|
||||
BRepLib::BuildCurves3d(edge);
|
||||
|
||||
wire_builder.Add(edge);
|
||||
wires.push_back (wire_builder.Wire());
|
||||
wire_builder = BRepBuilderAPI_MakeWire();
|
||||
return shared_from_this();
|
||||
}
|
||||
|
||||
auto NameVertex (string name)
|
||||
{
|
||||
if (!lastvertex.IsNull())
|
||||
@ -706,32 +689,6 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
|
||||
OCCGeometry::global_shape_properties.clear();
|
||||
OCCGeometry::global_shape_property_indices.Clear();
|
||||
});
|
||||
|
||||
struct SwigTypeInfo
|
||||
{
|
||||
const char* name; // SWIG's type name string
|
||||
// Other fields...
|
||||
};
|
||||
|
||||
struct SwigPyObject{
|
||||
PyObject_HEAD
|
||||
void *ptr;
|
||||
SwigTypeInfo* ty; // SWIG type information
|
||||
int own; // ownership flag
|
||||
};
|
||||
|
||||
m.def("From_PyOCC", [](py::object shape)
|
||||
{
|
||||
py::object py_this = shape.attr("this");
|
||||
PyObject* obj = py_this.ptr();
|
||||
SwigPyObject* swig_obj = reinterpret_cast<SwigPyObject*>(obj);
|
||||
if (!swig_obj->ptr || !swig_obj->ty || !swig_obj->ty->name) {
|
||||
throw std::runtime_error("SWIG object does not contain a valid pointer");
|
||||
}
|
||||
if(strcmp(swig_obj->ty->name, "_p_TopoDS_Shape") != 0)
|
||||
throw std::runtime_error("Does not contain TopoDS_Shape from pyocc!");
|
||||
return py::cast(static_cast<TopoDS_Shape*>(swig_obj->ptr));
|
||||
}, py::return_value_policy::reference, py::keep_alive<0,1>());
|
||||
|
||||
py::class_<TopoDS_Shape> (m, "TopoDS_Shape")
|
||||
.def("__str__", [] (const TopoDS_Shape & shape)
|
||||
@ -850,33 +807,10 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
|
||||
.def("WriteStep", [](const TopoDS_Shape & shape, string & filename)
|
||||
{ step_utils::WriteSTEP(shape, filename); }
|
||||
, py::arg("filename"), "export shape in STEP - format")
|
||||
.def("WriteBrep", [](const TopoDS_Shape & shape, const string& filename,
|
||||
bool withTriangles, bool withNormals,
|
||||
optional<int> version, bool binary)
|
||||
.def("WriteBrep", [](const TopoDS_Shape & shape, const string& filename)
|
||||
{
|
||||
if(binary)
|
||||
{
|
||||
#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
|
||||
BinTools_FormatVersion v = version ? BinTools_FormatVersion(*version) : BinTools_FormatVersion_CURRENT;
|
||||
BinTools::Write(shape, filename.c_str(), withTriangles, withNormals, v);
|
||||
# else // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
|
||||
throw Exception("Binary BREP export not supported in this version of OpenCascade");
|
||||
#endif // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
|
||||
}
|
||||
else
|
||||
{
|
||||
#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
|
||||
TopTools_FormatVersion v = version ? (TopTools_FormatVersion)(*version) : TopTools_FormatVersion_CURRENT;
|
||||
BRepTools::Write(shape, filename.c_str(), withTriangles, withNormals, v);
|
||||
#else // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
|
||||
BRepTools::Write(shape, filename.c_str());
|
||||
#endif // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
|
||||
}
|
||||
}, py::arg("filename"), py::arg("withTriangles")=true,
|
||||
py::arg("withNormals")=false,
|
||||
py::arg("version")=py::none(),
|
||||
py::arg("binary")=false,
|
||||
"export shape in BREP - format")
|
||||
BRepTools::Write(shape, filename.c_str());
|
||||
}, py::arg("filename"), "export shape in BREP - format")
|
||||
.def("bc", [](const TopoDS_Shape & shape, const string & name)
|
||||
{
|
||||
for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
|
||||
@ -1964,21 +1898,21 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
|
||||
})
|
||||
.def("Edge", [](Handle(Geom2d_Curve) curve) {
|
||||
// static Geom_Plane surf{gp_Ax3()}; // crashes in nbconvert ???
|
||||
static auto surf = Handle(Geom_Plane)(new Geom_Plane{gp_Ax3()});
|
||||
static auto surf = new Geom_Plane{gp_Ax3()};
|
||||
auto edge = BRepBuilderAPI_MakeEdge(curve, surf).Edge();
|
||||
BRepLib::BuildCurves3d(edge);
|
||||
return edge;
|
||||
})
|
||||
.def("Wire", [](Handle(Geom2d_Curve) curve) {
|
||||
// static Geom_Plane surf{gp_Ax3()}; // crashes in nbconvert ???
|
||||
static auto surf = Handle(Geom_Plane)(new Geom_Plane{gp_Ax3()});
|
||||
static auto surf = new Geom_Plane{gp_Ax3()};
|
||||
auto edge = BRepBuilderAPI_MakeEdge(curve, surf).Edge();
|
||||
BRepLib::BuildCurves3d(edge);
|
||||
return BRepBuilderAPI_MakeWire(edge).Wire();
|
||||
})
|
||||
.def("Face", [](Handle(Geom2d_Curve) curve) {
|
||||
// static Geom_Plane surf{gp_Ax3()}; // crashes in nbconvert ???
|
||||
static auto surf = Handle(Geom_Plane)(new Geom_Plane{gp_Ax3()});
|
||||
static auto surf = new Geom_Plane{gp_Ax3()};
|
||||
auto edge = BRepBuilderAPI_MakeEdge(curve, surf).Edge();
|
||||
BRepLib::BuildCurves3d(edge);
|
||||
auto wire = BRepBuilderAPI_MakeWire(edge).Wire();
|
||||
@ -2112,19 +2046,14 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
|
||||
return BRepOffsetAPI_MakePipe (spine, profile).Shape();
|
||||
}, py::arg("spine"), py::arg("profile"), py::arg("twist")=nullopt, py::arg("auxspine")=nullopt);
|
||||
|
||||
m.def("PipeShell", [] (const TopoDS_Wire & spine, variant<TopoDS_Shape, std::vector<TopoDS_Shape>> profile, std::optional<TopoDS_Wire> auxspine) {
|
||||
m.def("PipeShell", [] (const TopoDS_Wire & spine, const TopoDS_Shape & profile, const TopoDS_Wire & auxspine) {
|
||||
try
|
||||
{
|
||||
BRepOffsetAPI_MakePipeShell builder(spine);
|
||||
if(auxspine)
|
||||
builder.SetMode (*auxspine, Standard_True);
|
||||
if(std::holds_alternative<TopoDS_Shape>(profile))
|
||||
builder.Add (std::get<TopoDS_Shape>(profile));
|
||||
else
|
||||
{
|
||||
for(auto s : std::get<std::vector<TopoDS_Shape>>(profile))
|
||||
builder.Add(s);
|
||||
}
|
||||
builder.SetMode (auxspine, Standard_True);
|
||||
builder.Add (profile);
|
||||
// builder.Build();
|
||||
// builder.MakeSolid();
|
||||
return builder.Shape();
|
||||
}
|
||||
catch (Standard_Failure & e)
|
||||
@ -2133,13 +2062,13 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
|
||||
e.Print(errstr);
|
||||
throw NgException("cannot create PipeShell: "+errstr.str());
|
||||
}
|
||||
}, py::arg("spine"), py::arg("profile"), py::arg("auxspine")=nullopt);
|
||||
}, py::arg("spine"), py::arg("profile"), py::arg("auxspine"));
|
||||
|
||||
|
||||
// Handle(Geom2d_Ellipse) anEllipse1 = new Geom2d_Ellipse(anAx2d, aMajor, aMinor);
|
||||
m.def("Ellipse", [] (const gp_Ax2d & ax, double major, double minor) -> Handle(Geom2d_Curve)
|
||||
{
|
||||
return Handle(Geom2d_Ellipse) (GCE2d_MakeEllipse(ax, major, minor));
|
||||
return new Geom2d_Ellipse(ax, major, minor);
|
||||
}, py::arg("axes"), py::arg("major"), py::arg("minor"), "create 2d ellipse curve");
|
||||
|
||||
m.def("Segment", [](gp_Pnt2d p1, gp_Pnt2d p2) -> Handle(Geom2d_Curve) {
|
||||
@ -2707,8 +2636,6 @@ degen_tol : double
|
||||
.def("Circle", [](WorkPlane&wp, double x, double y, double r) {
|
||||
return wp.Circle(x,y,r); }, py::arg("h"), py::arg("v"), py::arg("r"), "draw circle with center (h,v) and radius 'r'")
|
||||
.def("Circle", [](WorkPlane&wp, double r) { return wp.Circle(r); }, py::arg("r"), "draw circle with center in current position")
|
||||
.def("Ellipse", [](WorkPlane& wp, double major, double minor)
|
||||
{ return wp.Ellipse(major, minor); }, py::arg("major"), py::arg("minor"), "draw ellipse with current position as center")
|
||||
.def("NameVertex", &WorkPlane::NameVertex, py::arg("name"), "name vertex at current position")
|
||||
.def("Offset", &WorkPlane::Offset, py::arg("d"), "replace current wire by offset curve of distance 'd'")
|
||||
.def("Reverse", &WorkPlane::Reverse, "revert orientation of current wire")
|
||||
|
||||
@ -57,13 +57,6 @@ namespace netgen
|
||||
return opengl_text_width;
|
||||
}
|
||||
|
||||
void MyOpenGLLines(FlatArray<Point<3>> points)
|
||||
{
|
||||
glBegin(GL_LINES);
|
||||
for (auto p : points)
|
||||
glVertex3dv(&p[0]);
|
||||
glEnd();
|
||||
}
|
||||
|
||||
// texture for color decoding
|
||||
// GLubyte * VisualScene :: colortexture = NULL;
|
||||
|
||||
@ -90,7 +90,6 @@ namespace netgen
|
||||
NGGUI_API extern void Set_OpenGLText_Callback ( void (*fun) (const char * text), int width );
|
||||
NGGUI_API extern VisualScene visual_scene_cross;
|
||||
NGGUI_API extern VisualScene *visual_scene;
|
||||
NGGUI_API extern void MyOpenGLLines (FlatArray<Point<3>> points);
|
||||
|
||||
|
||||
|
||||
|
||||
@ -841,7 +841,7 @@ namespace netgen
|
||||
if (vispar.clipping.enable && clipsolution == 2)
|
||||
{
|
||||
mesh->Mutex().unlock();
|
||||
mesh->BuildElementSearchTree(3);
|
||||
mesh->BuildElementSearchTree();
|
||||
mesh->Mutex().lock();
|
||||
}
|
||||
|
||||
@ -2667,8 +2667,6 @@ namespace netgen
|
||||
for (int i=first; i<next; i++)
|
||||
{
|
||||
double val;
|
||||
if(!VolumeElementActive(sol, *mesh, (*mesh)[ElementIndex(i)]))
|
||||
continue;
|
||||
bool considerElem = GetValue (sol, i, 0.333, 0.333, 0.333, comp, val);
|
||||
if (considerElem)
|
||||
{
|
||||
@ -2700,8 +2698,6 @@ namespace netgen
|
||||
// for (int i = 0; i < nse; i++)
|
||||
for (SurfaceElementIndex i : mesh->SurfaceElements().Range())
|
||||
{
|
||||
if(!SurfaceElementActive(sol, *mesh, (*mesh)[i]))
|
||||
continue;
|
||||
ELEMENT_TYPE type = (*mesh)[i].GetType();
|
||||
double val;
|
||||
bool considerElem = (type == QUAD)
|
||||
@ -4341,7 +4337,7 @@ namespace netgen
|
||||
}
|
||||
|
||||
double lami[3];
|
||||
int elnr = mesh->GetElementOfPoint (hp, lami,0,cindex,allowindex);
|
||||
int elnr = mesh->GetElementOfPoint (hp, lami,0,cindex,allowindex)-1;
|
||||
|
||||
if (elnr != -1)
|
||||
{
|
||||
@ -4727,7 +4723,7 @@ namespace netgen
|
||||
|
||||
|
||||
bool VisualSceneSolution ::
|
||||
SurfaceElementActive(const SolData *data, const Mesh & mesh, const Element2d & el) const
|
||||
SurfaceElementActive(const SolData *data, const Mesh & mesh, const Element2d & el)
|
||||
{
|
||||
if(data == nullptr) return true;
|
||||
bool is_active = true;
|
||||
@ -4753,7 +4749,7 @@ namespace netgen
|
||||
}
|
||||
|
||||
bool VisualSceneSolution ::
|
||||
VolumeElementActive(const SolData *data, const Mesh & mesh, const Element & el) const
|
||||
VolumeElementActive(const SolData *data, const Mesh & mesh, const Element & el)
|
||||
{
|
||||
bool is_active = true;
|
||||
if(data->draw_volumes)
|
||||
@ -4862,18 +4858,18 @@ namespace netgen
|
||||
if(sol.iscomplex && rcomponent != 0)
|
||||
{
|
||||
rcomponent = 2 * ((rcomponent-1)/2) + 1;
|
||||
GetValue(&sol, el3d, lami[0], lami[1], lami[2], rcomponent+1,
|
||||
GetValue(&sol, el3d-1, lami[0], lami[1], lami[2], rcomponent+1,
|
||||
imag);
|
||||
comp = (scalcomp-1)/2 + 1;
|
||||
}
|
||||
GetValue(&sol, el3d, lami[0], lami[1], lami[2], rcomponent, val);
|
||||
GetValue(&sol, el3d-1, lami[0], lami[1], lami[2], rcomponent, val);
|
||||
printScalValue(sol, comp, val, imag, sol.iscomplex && comp > 0);
|
||||
}
|
||||
if(vecfunction!=-1 && soldata[vecfunction]->draw_volume)
|
||||
{
|
||||
auto & sol = *soldata[vecfunction];
|
||||
ArrayMem<double, 10> values(sol.components);
|
||||
GetValues(&sol, el3d, lami[0], lami[1], lami[2], &values[0]);
|
||||
GetValues(&sol, el3d-1, lami[0], lami[1], lami[2], &values[0]);
|
||||
printVecValue(sol, values);
|
||||
}
|
||||
return;
|
||||
@ -4884,7 +4880,7 @@ namespace netgen
|
||||
double lami[3] = {0.0, 0.0, 0.0};
|
||||
// Check if unprojected Point is close to surface element (eps of 1e-3 due to z-Buffer accuracy)
|
||||
bool found_2del = false;
|
||||
if(selelement>0 && mesh->PointContainedIn2DElement(p, lami, selelement-1, false && fabs(lami[2])<1e-3))
|
||||
if(selelement>0 && mesh->PointContainedIn2DElement(p, lami, selelement, false && fabs(lami[2])<1e-3))
|
||||
{
|
||||
// Found it, use coordinates of point projected to surface element
|
||||
mesh->GetCurvedElements().CalcSurfaceTransformation({1.0-lami[0]-lami[1], lami[0]}, selelement-1, p);
|
||||
|
||||
@ -253,8 +253,8 @@ private:
|
||||
void DrawCone (const Point<3> & p1, const Point<3> & p2, double r);
|
||||
void DrawCylinder (const Point<3> & p1, const Point<3> & p2, double r);
|
||||
|
||||
bool SurfaceElementActive(const SolData *data, const Mesh & mesh, const Element2d & sei) const;
|
||||
bool VolumeElementActive(const SolData *data, const Mesh & mesh, const Element & ei) const;
|
||||
bool SurfaceElementActive(const SolData *data, const Mesh & mesh, const Element2d & sei);
|
||||
bool VolumeElementActive(const SolData *data, const Mesh & mesh, const Element & ei);
|
||||
|
||||
// Get Function Value, local coordinates lam1, lam2, lam3,
|
||||
bool GetValue (const SolData * data, ElementIndex elnr,
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user