ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 21:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

OCC Mesher Cleanup

Browse Source
This commit is contained in:
Matthias Hochsteger 2021-11-28 15:14:41 +00:00 committed by Joachim Schöberl
parent 2744d80aa3
commit 16b88e8e67
32 changed files with 2139 additions and 2043 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
cmake/SuperBuild.cmake
View File

@ -101,8 +101,8 @@ if(BUILD_OCC)
ExternalProject_Add(project_occ
DEPENDS project_freetype
URL https://github.com/Open-Cascade-SAS/OCCT/archive/refs/tags/V7_5_0.zip
URL_MD5 a24e6d3cf2d24bf9347d2d4aee9dd80a
URL https://github.com/Open-Cascade-SAS/OCCT/archive/refs/tags/V7_6_0.zip
URL_MD5 37519251c99cb3469ccfa82a9241d528
DOWNLOAD_DIR ${CMAKE_CURRENT_SOURCE_DIR}/external_dependencies
${SUBPROJECT_ARGS}
CMAKE_ARGS
@ -116,7 +116,7 @@ if(BUILD_OCC)
-DBUILD_MODULE_DataExchange:BOOL=ON
-DBUILD_MODULE_ApplicationFramework:BOOL=OFF
-DBUILD_MODULE_Draw:BOOL=OFF
-DUSE_FREETYPE=OFF
-DUSE_FREETYPE=ON
${SUBPROJECT_CMAKE_ARGS}
UPDATE_COMMAND ""
)

2
libsrc/gprim/adtree.hpp
View File

@ -817,6 +817,8 @@ public:
: BoxTree(box.PMin(), box.PMax())
{ }
double GetTolerance() { return tol; }
size_t GetNLeaves()
{
return n_leaves;

525
libsrc/meshing/basegeom.cpp
View File

@ -4,6 +4,28 @@
namespace netgen
{
struct PointTree
{
BoxTree<3> tree;
PointTree( Box<3> bb ) : tree(bb) {}
void Insert(Point<3> p, PointIndex n)
{
tree.Insert(p, p, n);
}
PointIndex Find(Point<3> p) const
{
ArrayMem<int, 1> points;
tree.GetIntersecting(p, p, points);
if(points.Size()==0)
throw Exception("cannot find mapped point");
return points[0];
}
double GetTolerance() { return tree.GetTolerance(); }
};
DLL_HEADER GeometryRegisterArray geometryregister;
//DLL_HEADER NgArray<GeometryRegister*> geometryregister;
@ -11,6 +33,38 @@ namespace netgen
GeometryRegister :: ~GeometryRegister()
{ ; }
bool GeometryShape :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const
{
throw Exception("GeometryShape::IsMappedShape not implemented for class " + Demangle(typeid(this).name()));
}
bool GeometryVertex :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const
{
const auto other_ptr = dynamic_cast<const GeometryVertex*>(&other_);
if(!other_ptr)
return false;
return Dist(trafo(GetPoint()), other_ptr->GetPoint()) < tol;
}
bool GeometryEdge :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const
{
const auto other_ptr = dynamic_cast<const GeometryEdge*>(&other_);
if(!other_ptr)
return false;
auto & e = *other_ptr;
if(tol < Dist(GetCenter(), e.GetCenter()))
return false;
auto &v0 = GetStartVertex();
auto &v1 = GetEndVertex();
auto &w0 = e.GetStartVertex();
auto &w1 = e.GetEndVertex();
return( (v0.IsMappedShape(w0, trafo, tol) && v1.IsMappedShape(w1, trafo, tol)) ||
(v0.IsMappedShape(w1, trafo, tol) && v1.IsMappedShape(w0, trafo, tol)) );
}
void GeometryFace :: RestrictHTrig(Mesh& mesh,
const PointGeomInfo& gi0,
const PointGeomInfo& gi1,
@ -103,6 +157,64 @@ namespace netgen
}
};
void NetgenGeometry :: ProcessIdentifications()
{
auto mirror_identifications = [&] ( auto & shapes )
{
for(auto i : Range(shapes))
{
auto &s = shapes[i];
s->nr = i;
for(auto & ident : s->identifications)
if(s.get() == ident.from)
ident.to->identifications.Append(ident);
}
};
mirror_identifications(vertices);
mirror_identifications(edges);
mirror_identifications(faces);
// todo: propagate identifications faces -> edges -> vertices
auto find_primary = [&] (auto & shapes)
{
for(auto &s : shapes)
{
s->primary = s.get();
s->primary_to_me = Transformation<3>{ Vec<3> {0,0,0} }; // init with identity
}
bool changed = true;
while(changed) {
changed = false;
for(auto &s : shapes)
{
auto current = s->primary;
for(auto & ident : current->identifications)
{
bool need_inverse = ident.from == s.get();
auto other = need_inverse ? ident.to : ident.from;
if(other->nr < current->nr)
{
auto trafo = ident.trafo;
if(need_inverse)
trafo = trafo.CalcInverse();
s->primary = other;
s->primary_to_me.Combine(trafo, s->primary_to_me);
changed = true;
}
}
}
}
};
find_primary(vertices);
find_primary(edges);
find_primary(faces);
}
void NetgenGeometry :: Analyse(Mesh& mesh,
const MeshingParameters& mparam) const
{
@ -241,60 +353,13 @@ namespace netgen
mesh.LoadLocalMeshSize(mparam.meshsizefilename);
}
void NetgenGeometry :: FindEdges(Mesh& mesh,
const MeshingParameters& mparam) const
void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params)
{
static Timer t1("MeshEdges"); RegionTimer regt(t1);
static Timer tdivide("Divide Edges");
static Timer tdivedgesections("Divide edge sections");
const char* savetask = multithread.task;
multithread.task = "Mesh Edges";
// create face descriptors and set bc names
mesh.SetNBCNames(faces.Size());
for(auto i : Range(faces.Size()))
{
mesh.SetBCName(i, faces[i]->GetName());
// todo find attached solids
FaceDescriptor fd(i+1, 1, 0, i+1);
fd.SetBCName(mesh.GetBCNamePtr(i));
mesh.AddFaceDescriptor(fd);
}
std::map<size_t, PointIndex> vert2meshpt;
for(auto i : Range(vertices))
{
const auto& vert = *vertices[i];
MeshPoint mp(vert.GetPoint());
vert2meshpt[vert.GetHash()] = mesh.AddPoint(mp);
}
size_t segnr = 0;
for(auto facenr : Range(faces.Size()))
{
const auto& face = *faces[facenr];
for(auto facebndnr : Range(face.GetNBoundaries()))
{
auto boundary = face.GetBoundary(facebndnr);
for(auto enr : Range(boundary))
{
multithread.percent = 100. * ((double(enr)/boundary.Size() + facebndnr)/face.GetNBoundaries() + facenr)/faces.Size();
const auto& oriented_edge = *boundary[enr];
auto edgenr = GetEdgeIndex(oriented_edge);
const auto& edge = edges[edgenr];
PointIndex startp, endp;
// throws if points are not found
startp = vert2meshpt.at(edge->GetStartVertex().GetHash());
endp = vert2meshpt.at(edge->GetEndVertex().GetHash());
// ignore collapsed edges
if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam())
continue;
Array<MeshPoint> mps;
Array<double> params;
static Timer tdivide("Divide Edges");
RegionTimer rt(tdivide);
// -------------------- DivideEdge -----------------
static constexpr size_t divide_edge_sections = 1000;
tdivide.Start();
double hvalue[divide_edge_sections+1];
hvalue[0] = 0;
@ -309,7 +374,7 @@ namespace netgen
}
int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
tdivedgesections.Stop();
mps.SetSize(nsubedges-1);
points.SetSize(nsubedges-1);
params.SetSize(nsubedges+1);
int i = 1;
@ -319,14 +384,14 @@ namespace netgen
if (hvalue[i1]/hvalue[divide_edge_sections]*nsubedges >= i)
{
params[i] = (double(i1)/divide_edge_sections);
mps[i-1] = MeshPoint(edge->GetPoint(params[i]));
points[i-1] = MeshPoint(edge->GetPoint(params[i]));
i++;
}
i1++;
if (i1 > divide_edge_sections)
{
nsubedges = i;
mps.SetSize(nsubedges-1);
points.SetSize(nsubedges-1);
params.SetSize(nsubedges+1);
cout << "divide edge: local h too small" << endl;
}
@ -339,32 +404,115 @@ namespace netgen
if(params[nsubedges] <= params[nsubedges-1])
{
cout << "CORRECTED" << endl;
mps.SetSize (nsubedges-2);
points.SetSize (nsubedges-2);
params.SetSize (nsubedges);
params[nsubedges-1] = 1.;
}
tdivide.Stop();
}
void NetgenGeometry :: FindEdges(Mesh& mesh,
const MeshingParameters& mparam) const
{
static Timer t1("MeshEdges"); RegionTimer regt(t1);
const char* savetask = multithread.task;
multithread.task = "Mesh Edges";
PointTree tree( bounding_box );
auto & identifications = mesh.GetIdentifications();
std::map<size_t, PointIndex> vert2meshpt;
for(auto & vert : vertices)
{
auto pi = mesh.AddPoint(vert->GetPoint());
tree.Insert(mesh[pi], pi);
vert2meshpt[vert->GetHash()] = pi;
mesh[pi].Singularity(vert->properties.hpref);
if(vert->properties.name)
{
Element0d el(pi, pi);
el.name = vert->properties.GetName();
mesh.SetCD3Name(pi, el.name);
mesh.pointelements.Append (el);
}
}
for(auto & vert : vertices)
for(auto & ident : vert->identifications)
identifications.Add(vert2meshpt[ident.from->GetHash()],
vert2meshpt[ident.to->GetHash()],
ident.name,
ident.type);
size_t segnr = 0;
auto nedges = edges.Size();
Array<Array<PointIndex>> all_pnums(nedges);
Array<Array<double>> all_params(nedges);
for (auto edgenr : Range(edges))
{
auto edge = edges[edgenr].get();
PointIndex startp, endp;
// throws if points are not found
startp = vert2meshpt.at(edge->GetStartVertex().GetHash());
endp = vert2meshpt.at(edge->GetEndVertex().GetHash());
// ignore collapsed edges
if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam())
continue;
// ----------- Add Points to mesh and create segments -----
Array<PointIndex> pnums(mps.Size() + 2);
auto & pnums = all_pnums[edgenr];
auto & params = all_params[edgenr];
Array<Point<3>> edge_points;
Array<double> edge_params;
if(edge->primary == edge)
{
DivideEdge(edge, mparam, mesh, edge_points, edge_params);
}
else
{
auto nr_primary = edge->primary->nr;
auto & pnums_primary = all_pnums[nr_primary];
auto & params_primary = all_params[nr_primary];
auto trafo = edge->primary_to_me;
auto np = pnums_primary.Size();
edge_points.SetSize(np-2);
edge_params.SetSize(np-2);
for(auto i : Range(np-2))
{
edge_points[i] = trafo(mesh[pnums_primary[i+1]]);
EdgePointGeomInfo gi;
edge->ProjectPoint(edge_points[i], &gi);
edge_params[i] = gi.dist;
}
// reverse entries if we have decreasing parameters
if(edge_params.Size()>2 && edge_params[0] > edge_params.Last())
for(auto i : Range((np-2)/2))
{
swap(edge_points[i], edge_points[np-3-i]);
swap(edge_params[i], edge_params[np-3-i]);
}
}
pnums.SetSize(edge_points.Size() + 2);
pnums[0] = startp;
pnums[mps.Size()+1] = endp;
pnums.Last() = endp;
double eps = bounding_box.Diam() * 1e-8;
params.SetSize(edge_points.Size()+2);
params[0] = 0.;
params.Last() = 1.;
for(auto i : Range(mps))
for(auto i : Range(edge_points))
{
bool exists = false;
for(auto pi : Range(mesh.Points()))
{
if((mesh[pi] - mps[i]).Length() < eps)
{
exists = true;
auto pi = mesh.AddPoint(edge_points[i]);
tree.Insert(mesh[pi], pi);
pnums[i+1] = pi;
break;
}
}
if(!exists)
pnums[i+1] = mesh.AddPoint(mps[i]);
params[i+1] = edge_params[i];
}
for(auto i : Range(pnums.Size()-1))
@ -373,28 +521,31 @@ namespace netgen
Segment seg;
seg[0] = pnums[i];
seg[1] = pnums[i+1];
seg.edgenr = segnr;
seg.edgenr = edgenr+1;
seg.si = edgenr+1;
seg.epgeominfo[0].dist = params[i];
seg.epgeominfo[1].dist = params[i+1];
seg.epgeominfo[0].edgenr = edgenr;
seg.epgeominfo[1].edgenr = edgenr;
seg.si = facenr+1;
seg.surfnr1 = facenr+1;
// TODO: implement functionality to transfer edge parameter t to face parameters u,v
for(auto j : Range(2))
face.CalcEdgePointGI(*edge, params[i+j],
seg.epgeominfo[j]);
if(!oriented_edge.OrientedLikeGlobal())
{
swap (seg[0], seg[1]);
swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
}
seg.singedge_left = edge->properties.hpref;
seg.singedge_right = edge->properties.hpref;
mesh.AddSegment(seg);
}
mesh.SetCD2Name(edgenr+1, edge->properties.GetName());
}
for (auto & edge : edges)
{
// identify points on edge
for(auto & ident : edge->identifications)
if(ident.from == edge.get())
{
auto & pnums = all_pnums[edge->nr];
// start and end vertex are already identified
for(auto pi : pnums.Range(1, pnums.Size()-1))
{
auto pi_other = tree.Find(ident.trafo(mesh[pi]));
identifications.Add(pi, pi_other, ident.name, ident.type);
}
}
}
@ -402,15 +553,8 @@ namespace netgen
multithread.task = savetask;
}
void NetgenGeometry :: MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) const
{
static Timer t1("Surface Meshing"); RegionTimer regt(t1);
const char* savetask = multithread.task;
multithread.task = "Mesh Surface";
Array<int, PointIndex> glob2loc(mesh.GetNP());
for(auto k : Range(faces))
bool NetgenGeometry :: MeshFace(Mesh& mesh, const MeshingParameters& mparam,
int k, FlatArray<int, PointIndex> glob2loc) const
{
multithread.percent = 100. * k/faces.Size();
const auto& face = *faces[k];
@ -420,9 +564,8 @@ namespace netgen
glob2loc = 0;
int cntp = 0;
for(auto& seg : mesh.LineSegments())
{
if(seg.si == k+1)
auto segments = face.GetBoundary(mesh);
for(auto& seg : segments)
{
for(auto j : Range(2))
{
@ -435,10 +578,7 @@ namespace netgen
}
}
}
}
for(auto & seg : mesh.LineSegments())
{
if(seg.si == k+1)
for(auto & seg : segments)
{
PointGeomInfo gi0, gi1;
gi0.trignum = gi1.trignum = k+1;
@ -450,7 +590,6 @@ namespace netgen
glob2loc[seg[1]],
gi0, gi1);
}
}
// TODO Set max area 2* area of face
@ -464,10 +603,170 @@ namespace netgen
{
mesh.SurfaceElements()[i].SetIndex(k+1);
}
return res != MESHING2_OK;
}
void NetgenGeometry :: MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) const
{
static Timer t1("Surface Meshing"); RegionTimer regt(t1);
const char* savetask = multithread.task;
multithread.task = "Mesh Surface";
size_t n_failed_faces = 0;
Array<int, PointIndex> glob2loc(mesh.GetNP());
for(auto k : Range(faces))
{
auto & face = *faces[k];
mesh.SetBCName(k, face.properties.GetName());
// todo find attached solids
FaceDescriptor fd(k+1, face.domin+1, face.domout+1, k+1);
fd.SetBCName(mesh.GetBCNamePtr(k));
mesh.AddFaceDescriptor(fd);
if(face.primary == &face)
{
if(MeshFace(mesh, mparam, k, glob2loc))
n_failed_faces++;
}
}
if(n_failed_faces)
{
cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl;
cout << "SURFACE MESHING ERROR OCCURRED IN " << n_failed_faces << " FACES:" << endl;
return;
}
if (mparam.perfstepsend >= MESHCONST_OPTSURFACE)
{
mesh.CalcSurfacesOfNode();
OptimizeSurface(mesh, mparam);
}
bool have_identifications = false;
for(auto & face : faces)
if(face->primary != face.get())
{
have_identifications = true;
MapSurfaceMesh(mesh, *face);
}
// identify points on faces
if(have_identifications)
{
mesh.CalcSurfacesOfNode();
BitArray is_identified_face(faces.Size());
is_identified_face = false;
for(auto & face : faces)
for(auto & ident : face->identifications)
{
is_identified_face.SetBit(ident.from->nr);
is_identified_face.SetBit(ident.to->nr);
}
PointTree tree( bounding_box );
Array<int, PointIndex> pi_to_face(mesh.GetNP());
pi_to_face = -1;
Array<SurfaceElementIndex> si_of_face;
Array<Array<PointIndex>> pi_of_face(faces.Size());
for(auto & face : faces)
if(is_identified_face[face->nr])
{
mesh.GetSurfaceElementsOfFace(face->nr+1, si_of_face);
for(auto si : si_of_face)
for(auto pi : mesh[si].PNums())
{
if(mesh[pi].Type() == SURFACEPOINT && pi_to_face[pi]==-1)
{
pi_to_face[pi] = face->nr;
tree.Insert(mesh[pi], pi);
pi_of_face[face->nr].Append(pi);
}
}
}
auto & mesh_ident = mesh.GetIdentifications();
for(auto & face : faces)
for(auto & ident : face->identifications)
{
if(ident.from == face.get())
for(auto pi : pi_of_face[face->nr])
{
auto pi_other = tree.Find(ident.trafo(mesh[pi]));
mesh_ident.Add(pi, pi_other, ident.name, ident.type);
}
}
}
mesh.CalcSurfacesOfNode();
multithread.task = savetask;
}
void NetgenGeometry :: MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const
{
static Timer timer("MapSurfaceMesh");
RegionTimer rt(timer);
const auto & src = dynamic_cast<const GeometryFace&>(*dst.primary);
auto trafo = dst.primary_to_me;
PrintMessage(2, "Map face ", src.nr+1, " -> ", dst.nr+1);
// point map from src to dst
Array<PointIndex, PointIndex> pmap(mesh.Points().Size());
pmap = PointIndex::INVALID;
// first map points on edges (mapped points alread in mesh, use search tree)
Array<bool, PointIndex> is_point_in_tree(mesh.Points().Size());
is_point_in_tree = false;
PointTree tree( bounding_box );
for (Segment & seg : src.GetBoundary(mesh))
for(auto i : Range(2))
{
auto pi = seg[i];
if(!is_point_in_tree[pi])
{
tree.Insert(trafo(mesh[pi]), pi);
is_point_in_tree[pi] = true;
}
}
for (Segment & seg : dst.GetBoundary(mesh))
for(auto i : Range(2))
{
auto pi = seg[i];
if(pmap[pi].IsValid())
continue;
pmap[tree.Find(mesh[pi])] = pi;
}
// now insert mapped surface elements
for(auto sei : mesh.SurfaceElements().Range())
{
auto sel = mesh[sei];
if(sel.GetIndex() != src.nr+1)
continue;
auto sel_new = sel;
sel_new.SetIndex(dst.nr+1);
for(auto i : Range(sel.PNums()))
{
auto pi = sel[i];
if(!pmap[pi].IsValid())
{
pmap[pi] = mesh.AddPoint(trafo(mesh[pi]), 1, SURFACEPOINT);
}
sel_new[i] = pmap[pi];
}
sel_new.Invert();
for(auto i : Range(sel.PNums()))
dst.CalcPointGeomInfo(mesh[sel_new[i]], sel_new.GeomInfo()[i]);
mesh.AddSurfaceElement(sel_new);
}
}
void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const
{
const auto savetask = multithread.task;
@ -505,6 +804,13 @@ namespace netgen
multithread.task = savetask;
}
void NetgenGeometry :: FinalizeMesh(Mesh& mesh) const
{
for (int i = 0; i < mesh.GetNDomains(); i++)
if (auto name = solids[i]->properties.name)
mesh.SetMaterial (i+1, *name);
}
shared_ptr<NetgenGeometry> GeometryRegisterArray :: LoadFromMeshFile (istream & ist) const
{
if (!ist.good())
@ -567,18 +873,17 @@ namespace netgen
if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
{
MeshSurface(*mesh, mparam);
mesh->CalcSurfacesOfNode();
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHSURFACE)
return 0;
if (mparam.perfstepsstart <= MESHCONST_OPTSURFACE)
OptimizeSurface(*mesh, mparam);
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE)
return 0;
if(dimension == 2)
{
FinalizeMesh(*mesh);
mesh->SetDimension(2);
return 0;
}
if(mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
{

93
libsrc/meshing/basegeom.hpp
View File

@ -12,26 +12,72 @@ struct Tcl_Interp;
namespace netgen
{
class GeometryVertex
struct ShapeProperties
{
public:
virtual ~GeometryVertex() {}
virtual Point<3> GetPoint() const = 0;
virtual size_t GetHash() const = 0;
optional<string> name;
optional<Vec<4>> col;
double maxh = 1e99;
double hpref = 0; // number of hp refinement levels (will be multiplied by factor later)
void Merge(const ShapeProperties & prop2)
{
if (prop2.name) name = prop2.name;
if (prop2.col) col = prop2.col;
maxh = min2(maxh, prop2.maxh);
hpref = max2(hpref, prop2.hpref);
}
string GetName() const { return name ? *name : "default"; }
Vec<4> GetColor() { return col ? *col : Vec<4>{0., 1., 0., 1.}; }
void DoArchive(Archive& ar)
{
ar & name & col & maxh & hpref;
}
};
class GeometryEdge
class GeometryShape;
struct ShapeIdentification
{
GeometryShape * from;
GeometryShape * to;
Transformation<3> trafo;
Identifications::ID_TYPE type;
string name = "";
};
class DLL_HEADER GeometryShape
{
public:
int nr = -1;
ShapeProperties properties;
Array<ShapeIdentification> identifications;
GeometryShape * primary;
Transformation<3> primary_to_me;
virtual ~GeometryShape() {}
virtual size_t GetHash() const = 0;
virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const;
};
class DLL_HEADER GeometryVertex : public GeometryShape
{
public:
virtual Point<3> GetPoint() const = 0;
virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override;
};
class DLL_HEADER GeometryEdge : public GeometryShape
{
public:
virtual ~GeometryEdge() {}
virtual const GeometryVertex& GetStartVertex() const = 0;
virtual const GeometryVertex& GetEndVertex() const = 0;
virtual double GetLength() const = 0;
virtual Point<3> GetCenter() const = 0;
virtual Point<3> GetPoint(double t) const = 0;
// Calculate parameter step respecting edges sag value
virtual double CalcStep(double t, double sag) const = 0;
virtual bool OrientedLikeGlobal() const = 0;
virtual size_t GetHash() const = 0;
virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const = 0;
virtual void PointBetween(const Point<3>& p1,
const Point<3>& p2,
@ -46,15 +92,17 @@ namespace netgen
ProjectPoint(newp, &newgi);
}
virtual Vec<3> GetTangent(double t) const = 0;
virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override;
};
class GeometryFace
class DLL_HEADER GeometryFace : public GeometryShape
{
public:
virtual ~GeometryFace() {}
int domin=-1, domout=-1;
virtual Point<3> GetCenter() const = 0;
virtual size_t GetNBoundaries() const = 0;
virtual Array<unique_ptr<GeometryEdge>> GetBoundary(size_t index) const = 0;
virtual string GetName() const { return "default"; }
virtual Array<Segment> GetBoundary(const Mesh& mesh) const = 0;
virtual PointGeomInfo Project(Point<3>& p) const = 0;
// Project point using geo info. Fast if point is close to
// parametrization in geo info.
@ -102,6 +150,9 @@ namespace netgen
int depth = 0, double h = 0.) const;
};
class DLL_HEADER GeometrySolid : public GeometryShape
{ };
class DLL_HEADER NetgenGeometry
{
unique_ptr<Refinement> ref;
@ -109,15 +160,25 @@ namespace netgen
Array<unique_ptr<GeometryVertex>> vertices;
Array<unique_ptr<GeometryEdge>> edges;
Array<unique_ptr<GeometryFace>> faces;
Array<unique_ptr<GeometrySolid>> solids;
Array<std::pair<Point<3>, double>> restricted_h;
Box<3> bounding_box;
int dimension = 3;
public:
NetgenGeometry()
{
ref = make_unique<Refinement>(*this);
}
virtual ~NetgenGeometry () { ; }
size_t GetNVertices() const { return vertices.Size(); }
size_t GetNEdges() const { return edges.Size(); }
size_t GetNFaces() const { return faces.Size(); }
const GeometryFace & GetFace(int i) const { return *faces[i]; }
virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
void RestrictH(const Point<3>& pnt, double maxh)
@ -134,13 +195,17 @@ namespace netgen
{ throw NgException("DoArchive not implemented for " + Demangle(typeid(*this).name())); }
virtual Mesh::GEOM_TYPE GetGeomType() const { return Mesh::NO_GEOM; }
virtual void ProcessIdentifications();
virtual void Analyse(Mesh& mesh,
const MeshingParameters& mparam) const;
virtual void FindEdges(Mesh& mesh, const MeshingParameters& mparam) const;
virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) const;
virtual bool MeshFace(Mesh& mesh, const MeshingParameters& mparam,
int nr, FlatArray<int, PointIndex> glob2loc) const;
virtual void MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const;
virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const;
virtual void FinalizeMesh(Mesh& mesh) const {}
virtual void FinalizeMesh(Mesh& mesh) const;
virtual PointGeomInfo ProjectPoint (int surfind, Point<3> & p) const
{

1
libsrc/meshing/meshtype.cpp
View File

@ -2687,6 +2687,7 @@ namespace netgen
identifiedpoints_nr.Set (tripl, 1);
if (identnr > maxidentnr) maxidentnr = identnr;
names.SetSize(maxidentnr);
if (identnr+1 > idpoints_table.Size())
idpoints_table.ChangeSize (identnr+1);

15
libsrc/meshing/meshtype.hpp
View File

@ -1497,6 +1497,7 @@ namespace netgen
/// number of identifications (or, actually used identifications ?)
int maxidentnr;
Array<string> names;
public:
///
@ -1511,7 +1512,12 @@ namespace netgen
identification nr identnr
*/
DLL_HEADER void Add (PointIndex pi1, PointIndex pi2, int identnr);
void Add (PointIndex pi1, PointIndex pi2, string name, ID_TYPE type)
{
auto nr = GetNr(name);
Add(pi1, pi2, nr);
SetType(nr, type);
}
int Get (PointIndex pi1, PointIndex pi2) const;
int GetSymmetric (PointIndex pi1, PointIndex pi2) const;
@ -1560,6 +1566,13 @@ namespace netgen
///
int GetMaxNr () const { return maxidentnr; }
int GetNr(string name)
{
if(!names.Contains(name))
names.Append(name);
return names.Pos(name)+1;
}
/// remove secondorder
void SetMaxPointNr (int maxpnum);

13
libsrc/meshing/python_mesh.cpp
View File

@ -117,6 +117,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
m.def("_PushStatus", [](string s) { PushStatus(MyStr(s)); });
m.def("_SetThreadPercentage", [](double percent) { SetThreadPercent(percent); });
py::enum_<Identifications::ID_TYPE>(m,"IdentificationType")
.value("UNDEFINED", Identifications::UNDEFINED)
.value("PERIODIC", Identifications::PERIODIC)
.value("CLOSESURFACES", Identifications::CLOSESURFACES)
.value("CLOSEEDGES", Identifications::CLOSEEDGES)
;
py::class_<NgMPI_Comm> (m, "MPI_Comm")
#ifdef NG_MPI4PY
.def(py::init([] (mpi4py_comm comm)
@ -944,19 +951,19 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def ("GetCD3Name", &Mesh::GetCD3Name)
.def ("SetCD3Name", &Mesh::SetCD3Name)
.def ("AddPointIdentification", [](Mesh & self, py::object pindex1, py::object pindex2, int identnr, int type)
.def ("AddPointIdentification", [](Mesh & self, py::object pindex1, py::object pindex2, int identnr, Identifications::ID_TYPE type)
{
if(py::extract<PointIndex>(pindex1).check() && py::extract<PointIndex>(pindex2).check())
{
self.GetIdentifications().Add (py::extract<PointIndex>(pindex1)(), py::extract<PointIndex>(pindex2)(), identnr);
self.GetIdentifications().SetType(identnr, Identifications::ID_TYPE(type)); // type = 2 ... periodic
self.GetIdentifications().SetType(identnr, type); // type = 2 ... periodic
}
},
//py::default_call_policies(),
py::arg("pid1"),
py::arg("pid2"),
py::arg("identnr"),
py::arg("type"))
py::arg("type")=Identifications::PERIODIC)
.def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries,
py::arg("face1"), py::arg("face2"), py::arg("mapping"), py::arg("point_tolerance") = -1.)
.def("GetNrIdentifications", [](Mesh& self)

12
libsrc/occ/CMakeLists.txt
View File

@ -1,9 +1,12 @@
if(USE_OCC)
add_definitions(-DNGINTERFACE_EXPORTS)
add_library(occ ${NG_LIB_TYPE}
Partition_Inter2d.cxx Partition_Inter3d.cxx
Partition_Loop.cxx Partition_Loop2d.cxx Partition_Loop3d.cxx Partition_Spliter.cxx
occgenmesh.cpp occgeom.cpp occmeshsurf.cpp python_occ.cpp
python_occ_basic.cpp python_occ_shapes.cpp
occ_face.cpp occ_edge.cpp occ_vertex.cpp occ_utils.cpp
)
if(USE_GUI)
add_library(occvis ${NG_LIB_TYPE} vsocc.cpp)
@ -14,11 +17,11 @@ target_link_libraries(occ PUBLIC ngcore PRIVATE "$<BUILD_INTERFACE:netgen_python
if(NOT WIN32)
target_link_libraries( occ PRIVATE ${OCC_LIBRARIES} )
if(USE_OCC AND APPLE)
if(APPLE)
# Link AppKit in case OCE was built as static libraries
find_library(AppKit AppKit)
target_link_libraries( occ PRIVATE ${AppKit} )
endif(USE_OCC AND APPLE)
endif(APPLE)
install( TARGETS occ ${NG_INSTALL_DIR})
if (USE_GUI)
target_link_libraries( occvis PUBLIC occ )
@ -27,6 +30,9 @@ if(NOT WIN32)
endif(NOT WIN32)
install(FILES
occgeom.hpp occmeshsurf.hpp vsocc.hpp
occgeom.hpp occmeshsurf.hpp vsocc.hpp occ_utils.hpp
occ_vertex.hpp occ_edge.hpp occ_face.hpp occ_solid.hpp
DESTINATION ${NG_INSTALL_DIR_INCLUDE}/occ COMPONENT netgen_devel
)
endif(USE_OCC)

84
libsrc/occ/occ_edge.cpp Normal file
View File

@ -0,0 +1,84 @@
#include <BRepGProp.hxx>
#include <BRep_Tool.hxx>
#include <GeomAPI_ProjectPointOnCurve.hxx>
#include "occ_edge.hpp"
#include "occgeom.hpp"
namespace netgen
{
OCCEdge::OCCEdge(TopoDS_Shape edge_)
: tedge(edge_.TShape()),
edge(TopoDS::Edge(edge_))
{
curve = BRep_Tool::Curve(edge, s0, s1);
BRepGProp::LinearProperties(edge, props);
auto verts = GetVertices(edge);
if(verts.size() != 2)
throw Exception("OCC edge does not have 2 vertices");
start = OCCVertex(verts[0]);
end = OCCVertex(verts[1]);
// swap start/end if necessary
double d00 = Dist(GetPoint(0), start.GetPoint());
double d01 = Dist(GetPoint(0), end.GetPoint());
if(d01 < d00)
swap(start, end);
}
const GeometryVertex& OCCEdge::GetStartVertex() const
{
return start;
}
const GeometryVertex& OCCEdge::GetEndVertex() const
{
return end;
}
double OCCEdge::GetLength() const
{
return props.Mass();
}
Point<3> OCCEdge::GetCenter() const
{
return occ2ng( props.CentreOfMass() );
}
Point<3> OCCEdge::GetPoint(double t) const
{
return occ2ng( curve->Value(s0+t*(s1-s0)) );
}
double OCCEdge::CalcStep(double t, double sag) const
{
throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
}
size_t OCCEdge::GetHash() const
{
return reinterpret_cast<size_t>(tedge.get());
}
void OCCEdge::ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const
{
auto pnt = ng2occ(p);
GeomAPI_ProjectPointOnCurve proj(pnt, curve);
pnt = proj.NearestPoint();
if(gi)
gi->dist = proj.LowerDistanceParameter();
p = occ2ng(pnt);
}
Vec<3> OCCEdge::GetTangent(double t) const
{
t = s0 + t*(s1-s0);
gp_Pnt p;
gp_Vec v;
curve->D1(t, p, v);
return occ2ng(v);
}
}

44
libsrc/occ/occ_edge.hpp Normal file
View File

@ -0,0 +1,44 @@
#ifndef FILE_OCC_EDGE_INCLUDED
#define FILE_OCC_EDGE_INCLUDED
#include <GProp_GProps.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Edge.hxx>
#include <Geom_Curve.hxx>
#include <BRep_TEdge.hxx>
#include "occ_vertex.hpp"
#include "meshing.hpp"
namespace netgen
{
class OCCEdge : public GeometryEdge
{
T_Shape tedge;
TopoDS_Edge edge;
Handle(Geom_Curve) curve;
double s0, s1;
GProp_GProps props;
OCCVertex start;
OCCVertex end;
public:
OCCEdge(TopoDS_Shape edge_);
auto Shape() const { return edge; }
T_Shape TShape() const { return tedge; }
const GeometryVertex& GetStartVertex() const override;
const GeometryVertex& GetEndVertex() const override;
double GetLength() const override;
Point<3> GetCenter() const override;
Point<3> GetPoint(double t) const override;
double CalcStep(double t, double sag) const override;
size_t GetHash() const override;
void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const override;
Vec<3> GetTangent(double t) const override;
};
}
#endif // FILE_OCCEDGE_INCLUDED

254
libsrc/occ/occ_face.cpp Normal file
View File

@ -0,0 +1,254 @@
#include <BRepGProp.hxx>
#include <BRep_Tool.hxx>
#include <GeomAPI_ProjectPointOnCurve.hxx>
#include "occ_edge.hpp"
#include "occ_face.hpp"
#include "occgeom.hpp"
namespace netgen
{
OCCFace::OCCFace(TopoDS_Shape dshape)
: tface(dshape.TShape()),
face(TopoDS::Face(dshape))
{
BRepGProp::LinearProperties(face, props);
bbox = ::netgen::GetBoundingBox(face);
surface = BRep_Tool::Surface(face);
shape_analysis = new ShapeAnalysis_Surface( surface );
tolerance = BRep_Tool::Tolerance( face );
}
size_t OCCFace::GetNBoundaries() const
{
return 0;
}
size_t OCCFace::GetHash() const
{
return reinterpret_cast<size_t>(tface.get());
}
Point<3> OCCFace::GetCenter() const
{
return occ2ng( props.CentreOfMass() );
}
Array<Segment> OCCFace::GetBoundary(const Mesh& mesh) const
{
auto & geom = dynamic_cast<OCCGeometry&>(*mesh.GetGeometry());
auto n_edges = geom.edge_map.size();
constexpr int UNUSED = 0;
constexpr int FORWARD = 1;
constexpr int REVERSED = 2;
constexpr int BOTH = 3;
Array<int> edge_orientation(n_edges);
edge_orientation = UNUSED;
Array<Handle(Geom2d_Curve)> curve_on_face[BOTH];
curve_on_face[FORWARD].SetSize(n_edges);
curve_on_face[REVERSED].SetSize(n_edges);
Array<TopoDS_Edge> edge_on_face[BOTH];
edge_on_face[FORWARD].SetSize(n_edges);
edge_on_face[REVERSED].SetSize(n_edges);
for(auto edge_ : GetEdges(face))
{
auto edge = TopoDS::Edge(edge_);
if(geom.edge_map.count(edge.TShape())==0)
continue;
auto edgenr = geom.edge_map[edge.TShape()];
auto & orientation = edge_orientation[edgenr];
double s0, s1;
auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
if(edge.Orientation() == TopAbs_FORWARD)
{
curve_on_face[FORWARD][edgenr] = cof;
orientation += FORWARD;
edge_on_face[FORWARD][edgenr] = edge;
}
if(edge.Orientation() == TopAbs_REVERSED)
{
curve_on_face[REVERSED][edgenr] = cof;
orientation += REVERSED;
edge_on_face[REVERSED][edgenr] = edge;
}
if(orientation > BOTH)
throw Exception("have edge more than twice in face " + ToString(nr) + " " + properties.GetName() + ", orientation: " + ToString(orientation));
}
Array<Segment> boundary;
for (auto seg : mesh.LineSegments())
{
auto edgenr = seg.epgeominfo[0].edgenr;
auto orientation = edge_orientation[edgenr];
if(orientation == UNUSED)
continue;
for(const auto ORIENTATION : {FORWARD, REVERSED})
{
if((orientation & ORIENTATION) == 0)
continue;
// auto cof = curve_on_face[ORIENTATION][edgenr];
auto edge = edge_on_face[ORIENTATION][edgenr];
OCCEdge gedge(edge);
double s0, s1;
auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
for(auto i : Range(2))
{
Point<3> p = mesh[seg[i]];
gedge.ProjectPoint(p, &seg.epgeominfo[i]);
}
double s[2] = { seg.epgeominfo[0].dist, seg.epgeominfo[1].dist };
// fixes normal-vector roundoff problem when endpoint is cone-tip
double delta = s[1]-s[0];
s[0] += 1e-10*delta;
s[1] -= 1e-10*delta;
for(auto i : Range(2))
{
auto uv = cof->Value(s[i]);
seg.epgeominfo[i].u = uv.X();
seg.epgeominfo[i].v = uv.Y();
}
if(ORIENTATION == REVERSED)
{
swap(seg[0], seg[1]);
swap(seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
swap(seg.epgeominfo[0].u, seg.epgeominfo[1].u);
swap(seg.epgeominfo[0].v, seg.epgeominfo[1].v);
}
boundary.Append(seg);
}
}
return boundary;
}
PointGeomInfo OCCFace::Project(Point<3>& p) const
{
auto suval = shape_analysis->ValueOfUV(ng2occ(p), tolerance);
double u,v;
suval.Coord(u, v);
p = occ2ng(surface->Value( u, v ));
PointGeomInfo gi;
gi.trignum = nr+1;
gi.u = u;
gi.v = v;
return gi;
}
bool OCCFace::ProjectPointGI(Point<3>& p_, PointGeomInfo& gi) const
{
double u = gi.u;
double v = gi.v;
auto p = ng2occ(p_);
auto x = surface->Value (u,v);
if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
gp_Vec du, dv;
surface->D1(u,v,x,du,dv);
int count = 0;
gp_Pnt xold;
gp_Vec n;
double det, lambda, mu;
do {
count++;
n = du^dv;
det = Det3 (n.X(), du.X(), dv.X(),
n.Y(), du.Y(), dv.Y(),
n.Z(), du.Z(), dv.Z());
if (det < 1e-15) return false;
lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
n.Y(), p.Y()-x.Y(), dv.Y(),
n.Z(), p.Z()-x.Z(), dv.Z())/det;
mu = Det3 (n.X(), du.X(), p.X()-x.X(),
n.Y(), du.Y(), p.Y()-x.Y(),
n.Z(), du.Z(), p.Z()-x.Z())/det;
u += lambda;
v += mu;
xold = x;
surface->D1(u,v,x,du,dv);
} while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
// (*testout) << "FastProject count: " << count << endl;
if (count == 50) return false;
p_ = occ2ng(x);
return true;
}
Point<3> OCCFace::GetPoint(const PointGeomInfo& gi) const
{
return occ2ng(surface->Value( gi.u, gi.v ));
}
void OCCFace::CalcEdgePointGI(const GeometryEdge& edge,
double t,
EdgePointGeomInfo& egi) const
{
throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
}
Box<3> OCCFace::GetBoundingBox() const
{
return bbox;
}
double OCCFace::GetCurvature(const PointGeomInfo& gi) const
{
throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
}
void OCCFace::RestrictH(Mesh& mesh, const MeshingParameters& mparam) const
{
throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
}
Vec<3> OCCFace::GetNormal(const Point<3>& p, const PointGeomInfo* gi) const
{
PointGeomInfo gi_;
if(gi==nullptr)
{
auto p_ = p;
gi_ = Project(p_);
gi = &gi_;
}
gp_Pnt pnt;
gp_Vec du, dv;
surface->D1(gi->u,gi->v,pnt,du,dv);
auto n = Cross (occ2ng(du), occ2ng(dv));
n.Normalize();
if (face.Orientation() == TopAbs_REVERSED)
n *= -1;
return n;
}
}

49
libsrc/occ/occ_face.hpp Normal file
View File

@ -0,0 +1,49 @@
#ifndef FILE_OCC_FACE_INCLUDED
#define FILE_OCC_FACE_INCLUDED
#include <GProp_GProps.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Face.hxx>
#include <ShapeAnalysis_Surface.hxx>
#include "occ_vertex.hpp"
#include "meshing.hpp"
namespace netgen
{
class OCCFace : public GeometryFace
{
T_Shape tface;
TopoDS_Face face;
GProp_GProps props;
Box<3> bbox;
Handle( Geom_Surface ) surface;
Handle( ShapeAnalysis_Surface ) shape_analysis;
double tolerance;
public:
OCCFace(TopoDS_Shape dshape);
T_Shape TShape() { return tface; }
size_t GetHash() const override;
Point<3> GetCenter() const override;
virtual size_t GetNBoundaries() const override;
virtual Array<Segment> GetBoundary(const Mesh& mesh) const override;
virtual PointGeomInfo Project(Point<3>& p) const override;
virtual bool ProjectPointGI(Point<3>& p, PointGeomInfo& gi) const override;
virtual Point<3> GetPoint(const PointGeomInfo& gi) const override;
virtual void CalcEdgePointGI(const GeometryEdge& edge,
double t,
EdgePointGeomInfo& egi) const override;
virtual Box<3> GetBoundingBox() const override;
virtual double GetCurvature(const PointGeomInfo& gi) const override;
virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const override;
virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo* gi = nullptr) const override;
};
}
#endif // FILE_OCC_FACE_INCLUDED

27
libsrc/occ/occ_solid.hpp Normal file
View File

@ -0,0 +1,27 @@
#ifndef FILE_OCC_SOLID_INCLUDED
#define FILE_OCC_SOLID_INCLUDED
#include <TopoDS.hxx>
#include <TopoDS_Solid.hxx>
#include "meshing.hpp"
namespace netgen
{
class OCCSolid : public GeometrySolid
{
T_Shape tsolid;
TopoDS_Solid solid;
public:
OCCSolid(TopoDS_Shape dshape)
: tsolid(dshape.TShape()),
solid(TopoDS::Solid(dshape))
{ }
T_Shape TShape() { return tsolid; }
size_t GetHash() const override { return reinterpret_cast<size_t>(tsolid.get()); }
};
}
#endif // FILE_OCC_SOLID_INCLUDED

24
libsrc/occ/occ_utils.cpp Normal file
View File

@ -0,0 +1,24 @@
#include <Bnd_Box.hxx>
#include <BRepBndLib.hxx>
#include <BRep_TVertex.hxx>
#include "occ_utils.hpp"
namespace netgen
{
Point<3> occ2ng (Handle(TopoDS_TShape) shape)
{
return occ2ng( Handle(BRep_TVertex)::DownCast(shape)->Pnt() );
}
Box<3> GetBoundingBox( const TopoDS_Shape & shape )
{
Bnd_Box bb;
#if OCC_VERSION_HEX < 0x070000
BRepBndLib::Add (shape, bb);
#else
BRepBndLib::Add (shape, bb, true);
#endif
return {occ2ng(bb.CornerMin()), occ2ng(bb.CornerMax())};
}
}

267
libsrc/occ/occ_utils.hpp Normal file
View File

@ -0,0 +1,267 @@
#ifndef FILE_OCC_UTILS_INCLUDED
#define FILE_OCC_UTILS_INCLUDED
#include <BRepGProp.hxx>
#include <BRep_Tool.hxx>
#include <GProp_GProps.hxx>
#include <TopExp_Explorer.hxx>
#include <TopTools_IndexedMapOfShape.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Vertex.hxx>
#include "meshing.hpp"
#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
#define OCC_HAVE_DUMP_JSON
#endif
namespace netgen
{
typedef Handle(TopoDS_TShape) T_Shape;
inline Point<3> occ2ng (const gp_Pnt & p)
{
return Point<3> (p.X(), p.Y(), p.Z());
}
inline Point<2> occ2ng (const gp_Pnt2d & p)
{
return Point<2> (p.X(), p.Y());
}
inline Vec<3> occ2ng (const gp_Vec & v)
{
return Vec<3> (v.X(), v.Y(), v.Z());
}
DLL_HEADER Point<3> occ2ng (T_Shape shape);
inline Point<3> occ2ng (const TopoDS_Shape & s)
{
return occ2ng(s.TShape());
}
inline Point<3> occ2ng (const TopoDS_Vertex & v)
{
return occ2ng (BRep_Tool::Pnt (v));
}
inline gp_Pnt ng2occ (const Point<3> & p)
{
return gp_Pnt(p(0), p(1), p(2));
}
DLL_HEADER Box<3> GetBoundingBox( const TopoDS_Shape & shape );
class OCCIdentification
{
public:
T_Shape from;
T_Shape to;
Transformation<3> trafo;
string name;
Identifications::ID_TYPE type;
bool opposite_direction;
};
class MyExplorer
{
class Iterator
{
TopExp_Explorer exp;
public:
Iterator (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid)
: exp(ashape, atoFind, atoAvoid) { }
auto operator*() { return exp.Current(); }
Iterator & operator++() { exp.Next(); return *this; }
bool operator!= (nullptr_t nu) { return exp.More(); }
};
public:
TopoDS_Shape shape;
TopAbs_ShapeEnum toFind;
TopAbs_ShapeEnum toAvoid;
MyExplorer (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid = TopAbs_SHAPE)
: shape(ashape), toFind(atoFind), toAvoid(atoAvoid) { ; }
Iterator begin() { return Iterator(shape, toFind, toAvoid); }
auto end() { return nullptr; }
};
inline auto Explore (TopoDS_Shape shape, TopAbs_ShapeEnum toFind, TopAbs_ShapeEnum toAvoid = TopAbs_SHAPE)
{
return MyExplorer (shape, toFind, toAvoid);
}
class IndexMapIterator
{
class Iterator
{
const TopTools_IndexedMapOfShape & indmap;
int i;
public:
Iterator (const TopTools_IndexedMapOfShape & aindmap, int ai)
: indmap(aindmap), i(ai) { ; }
auto operator*() { return tuple(i, indmap(i)); }
Iterator & operator++() { i++; return *this; }
bool operator!= (const Iterator & i2) { return i != i2.i; }
};
public:
const TopTools_IndexedMapOfShape & indmap;
IndexMapIterator (const TopTools_IndexedMapOfShape & aindmap) : indmap(aindmap) { }
Iterator begin() { return Iterator(indmap, 1); }
Iterator end() { return Iterator(indmap, indmap.Extent()+1); }
};
inline auto Enumerate (const TopTools_IndexedMapOfShape & indmap)
{
return IndexMapIterator(indmap);
}
struct ShapeLess
{
bool operator() (const TopoDS_Shape& s1, const TopoDS_Shape& s2) const
{
return s1.TShape() < s2.TShape();
}
};
class ListOfShapes : public std::vector<TopoDS_Shape>
{
public:
DLL_HEADER TopoDS_Shape Max(gp_Vec dir);
DLL_HEADER TopoDS_Shape Nearest(gp_Pnt pnt);
DLL_HEADER ListOfShapes SubShapes(TopAbs_ShapeEnum type) const;
ListOfShapes Solids() const
{
return SubShapes(TopAbs_SOLID);
}
ListOfShapes Faces() const
{
return SubShapes(TopAbs_FACE);
}
ListOfShapes Edges() const
{
return SubShapes(TopAbs_EDGE);
}
ListOfShapes Vertices() const
{
return SubShapes(TopAbs_VERTEX);
}
ListOfShapes operator*(const ListOfShapes& other) const
{
ListOfShapes common;
for(const auto& shape : (*this))
for(const auto& shape_o : other)
if(shape.IsSame(shape_o))
common.push_back(shape);
return common;
}
};
inline ListOfShapes GetSolids(const TopoDS_Shape & shape)
{
ListOfShapes sub;
for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
sub.push_back(e.Current());
return sub;
}
inline ListOfShapes GetFaces(const TopoDS_Shape & shape)
{
ListOfShapes sub;
for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
sub.push_back(e.Current());
return sub;
}
inline ListOfShapes GetEdges(const TopoDS_Shape & shape)
{
ListOfShapes sub;
for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next())
sub.push_back(e.Current());
return sub;
}
inline ListOfShapes GetVertices(const TopoDS_Shape & shape)
{
ListOfShapes sub;
for (TopExp_Explorer e(shape, TopAbs_VERTEX); e.More(); e.Next())
sub.push_back(e.Current());
return sub;
}
class DirectionalInterval
{
public:
gp_Vec dir;
double minval = -1e99;
double maxval = 1e99;
bool openmin = false, openmax = false;
DirectionalInterval (gp_Vec adir) : dir(adir) { ; }
DirectionalInterval (const DirectionalInterval & i2)
: dir(i2.dir), minval(i2.minval), maxval(i2.maxval) { ; }
DirectionalInterval operator< (double val) const
{
DirectionalInterval i2 = *this;
i2.maxval = val;
return i2;
}
DirectionalInterval operator> (double val) const
{
DirectionalInterval i2 = *this;
i2.minval = val;
return i2;
}
DirectionalInterval Intersect (const DirectionalInterval & i2)
{
DirectionalInterval res = *this;
res.minval = max(res.minval, i2.minval);
res.maxval = min(res.maxval, i2.maxval);
return res;
}
bool Contains (gp_Pnt p, double eps = 1e-8)
{
// cout << "Contains point " << p.X() << "," << p.Y() << "," << p.Z() << " ? " << endl;
double val = dir.X()*p.X() + dir.Y()*p.Y() + dir.Z() * p.Z();
// cout << "minval = " << minval << ", val = " << val << " maxval = " << maxval << endl;
if (openmin) {
if (val < minval+eps) return false;
} else {
if (val < minval-eps) return false;
}
if (openmax) {
if (val > maxval-eps) return false;
} else {
if (val > maxval+eps) return false;
}
return true;
}
};
inline gp_Pnt Center (TopoDS_Shape shape)
{
GProp_GProps props;
switch (shape.ShapeType())
{
case TopAbs_FACE:
BRepGProp::SurfaceProperties (shape, props); break;
default:
BRepGProp::LinearProperties(shape, props);
}
return props.CentreOfMass();
}
}
#endif // FILE_OCC_UTILS_INCLUDED

25
libsrc/occ/occ_vertex.cpp Normal file
View File

@ -0,0 +1,25 @@
#include <BRepGProp.hxx>
#include <BRep_Tool.hxx>
#include "occ_vertex.hpp"
namespace netgen
{
OCCVertex::OCCVertex( TopoDS_Shape s )
: vertex(TopoDS::Vertex(s)),
tvertex(s.TShape())
{
p = occ2ng(vertex);
}
Point<3> OCCVertex::GetPoint() const
{
return p;
}
size_t OCCVertex::GetHash() const
{
return reinterpret_cast<size_t>(tvertex.get());
}
}

28
libsrc/occ/occ_vertex.hpp Normal file
View File

@ -0,0 +1,28 @@
#ifndef FILE_OCC_VERTEX_INCLUDED
#define FILE_OCC_VERTEX_INCLUDED
#include <TopoDS.hxx>
#include <BRep_TVertex.hxx>
#include "meshing.hpp"
#include "occ_utils.hpp"
namespace netgen
{
class OCCVertex : public GeometryVertex
{
TopoDS_Vertex vertex;
T_Shape tvertex;
Point<3> p;
public:
OCCVertex( ) = default;
OCCVertex( TopoDS_Shape s );
~OCCVertex() {}
Point<3> GetPoint() const override;
size_t GetHash() const override;
T_Shape TShape() { return tvertex; }
};
}
#endif // FILE_OCC_VERTEX_INCLUDED

730
libsrc/occ/occgenmesh.cpp
View File

@ -1,14 +1,27 @@
#ifdef OCCGEOMETRY
#include <mystdlib.h>
#include <occgeom.hpp>
#include <meshing.hpp>
#include "occgeom.hpp"
#include "occmeshsurf.hpp"
#include <BRepAdaptor_Curve.hxx>
#include <BRepGProp.hxx>
#include <BRepLProp_CLProps.hxx>
#include <BRepLProp_SLProps.hxx>
#include <BRepMesh_IncrementalMesh.hxx>
#include <BRepTools.hxx>
#include <GProp_GProps.hxx>
#include <Quantity_Color.hxx>
#include <ShapeAnalysis.hxx>
#include <TopExp.hxx>
#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
#include <TopoDS_Edge.hxx>
namespace netgen
{
#include "occmeshsurf.hpp"
#define TCL_OK 0
#define TCL_ERROR 1
@ -218,536 +231,27 @@ namespace netgen
}
}
void DivideEdge (TopoDS_Edge & edge, NgArray<MeshPoint> & ps,
Array<double> & params, Mesh & mesh,
const MeshingParameters & mparam)
{
double s0, s1;
int nsubedges = 1;
gp_Pnt pnt, oldpnt;
double svalue[DIVIDEEDGESECTIONS];
GProp_GProps system;
BRepGProp::LinearProperties(edge, system);
double L = system.Mass();
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
double hvalue[DIVIDEEDGESECTIONS+1];
hvalue[0] = 0;
pnt = c->Value(s0);
int tmpVal = (int)(DIVIDEEDGESECTIONS);
for (int i = 1; i <= tmpVal; i++)
{
oldpnt = pnt;
pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));
hvalue[i] = hvalue[i-1] +
1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
pnt.Distance(oldpnt);
//(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))
// << " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;
}
// nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS]));
nsubedges = max (1, int(floor(hvalue[DIVIDEEDGESECTIONS]+0.5)));
ps.SetSize(nsubedges-1);
params.SetSize(nsubedges+1);
int i = 1;
int i1 = 0;
do
{
if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i)
{
params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0);
pnt = c->Value(params[i]);
ps[i-1] = MeshPoint (Point3d(pnt.X(), pnt.Y(), pnt.Z()));
i++;
}
i1++;
if (i1 > DIVIDEEDGESECTIONS)
{
nsubedges = i;
ps.SetSize(nsubedges-1);
params.SetSize(nsubedges+1);
cout << "divide edge: local h too small" << endl;
}
} while (i < nsubedges);
params[0] = s0;
params[nsubedges] = s1;
if (params[nsubedges] <= params[nsubedges-1])
{
cout << "CORRECTED" << endl;
ps.SetSize (nsubedges-2);
params.SetSize (nsubedges);
params[nsubedges] = s1;
}
}
void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam)
{
static Timer t("OCCFindEdges"); RegionTimer r(t);
static Timer tsearch("OCCFindEdges - search point");
const char * savetask = multithread.task;
multithread.task = "Edge meshing";
(*testout) << "edge meshing" << endl;
int nvertices = geom.vmap.Extent();
int nedges = geom.emap.Extent();
Array<Array<PointIndex>> alledgepnums(nedges);
Array<Array<double>> alledgeparams(nedges);
(*testout) << "nvertices = " << nvertices << endl;
(*testout) << "nedges = " << nedges << endl;
double eps = 1e-6 * geom.GetBoundingBox().Diam();
tsearch.Start();
for (auto [i,vshape] : Enumerate(geom.vmap))
{
TopoDS_Vertex vertex = TopoDS::Vertex(vshape);
gp_Pnt pnt = BRep_Tool::Pnt (vertex);
mesh.AddPoint (occ2ng(pnt));
double hpref = OCCGeometry::global_shape_properties[vertex.TShape()].hpref;
mesh.Points().Last().Singularity(hpref);
double maxh = OCCGeometry::global_shape_properties[vertex.TShape()].maxh;
mesh.RestrictLocalH (occ2ng(pnt), maxh);
}
tsearch.Stop();
(*testout) << "different vertices = " << mesh.GetNP() << endl;
// int first_ep = mesh.GetNP()+1;
// PointIndex first_ep = mesh.Points().End();
PointIndex first_ep = *mesh.Points().Range().end();
auto vertexrange = mesh.Points().Range();
NgArray<int> face2solid[2];
for (int i = 0; i < 2; i++)
{
face2solid[i].SetSize (geom.fmap.Extent());
face2solid[i] = 0;
}
/*
int solidnr = 0;
for (TopExp_Explorer exp0(geom.shape, TopAbs_SOLID); exp0.More(); exp0.Next())
{
solidnr++;
for (TopExp_Explorer exp1(exp0.Current(), TopAbs_FACE); exp1.More(); exp1.Next())
{
TopoDS_Face face = TopoDS::Face(exp1.Current());
int facenr = geom.fmap.FindIndex(face);
if(facenr < 1) continue;
if (face2solid[0][facenr-1] == 0)
face2solid[0][facenr-1] = solidnr;
else
face2solid[1][facenr-1] = solidnr;
}
}
*/
int solidnr = 0;
for (auto solid : Explore(geom.shape, TopAbs_SOLID))
{
solidnr++;
for (auto face : Explore (solid, TopAbs_FACE))
if (geom.fmap.Contains(face))
{
int facenr = geom.fmap.FindIndex(face);
if (face2solid[0][facenr-1] == 0)
face2solid[0][facenr-1] = solidnr;
else
face2solid[1][facenr-1] = solidnr;
}
}
/*
int total = 0;
for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++)
for (TopExp_Explorer exp2(geom.fmap(i3), TopAbs_WIRE); exp2.More(); exp2.Next())
for (TopExp_Explorer exp3(exp2.Current(), TopAbs_EDGE); exp3.More(); exp3.Next())
total++;
*/
int total = 0;
for (auto [i3, face] : Enumerate(geom.fmap))
for (auto wire : Explore(face, TopAbs_WIRE))
for (auto edge : Explore(wire, TopAbs_EDGE))
total++;
int facenr = 0;
// int edgenr = mesh.GetNSeg();
(*testout) << "faces = " << geom.fmap.Extent() << endl;
int curr = 0;
/*
for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++)
{
TopoDS_Face face = TopoDS::Face(geom.fmap(i3));
*/
for (auto [i3,faceshape] : Enumerate(geom.fmap))
{
TopoDS_Face face = TopoDS::Face(faceshape);
facenr = geom.fmap.FindIndex (face); // sollte doch immer == i3 sein ??? JS
if (facenr != i3)
cout << "info: facenr != i3, no problem, but please report to developers" << endl;
int solidnr0 = face2solid[0][i3-1];
int solidnr1 = face2solid[1][i3-1];
/* auskommentiert am 3.3.05 von robert
for (exp2.Init (geom.somap(solidnr0), TopAbs_FACE); exp2.More(); exp2.Next())
{
TopoDS_Face face2 = TopoDS::Face(exp2.Current());
if (geom.fmap.FindIndex(face2) == facenr)
{
// if (face.Orientation() != face2.Orientation()) swap (solidnr0, solidnr1);
}
}
*/
mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0));
Vec<4> col = OCCGeometry::global_shape_properties[face.TShape()].col.value_or(Vec<4>(0,1,0,1));
mesh.GetFaceDescriptor(facenr).SetSurfColour(col);
if(auto & opt_name = geom.fprops[facenr-1]->name)
mesh.GetFaceDescriptor(facenr).SetBCName(*opt_name);
else
mesh.GetFaceDescriptor(facenr).SetBCName("bc_"+ToString(facenr));
mesh.GetFaceDescriptor(facenr).SetBCProperty(facenr);
// ACHTUNG! STIMMT NICHT ALLGEMEIN (RG)
// kA was RG damit meinte
Handle(Geom_Surface) occface = BRep_Tool::Surface(face);
/*
for (TopExp_Explorer exp2 (face, TopAbs_WIRE); exp2.More(); exp2.Next())
{
TopoDS_Shape wire = exp2.Current();
*/
for (auto wire : MyExplorer (face, TopAbs_WIRE))
{
// for (TopExp_Explorer exp3 (wire, TopAbs_EDGE); exp3.More(); exp3.Next())
for (auto edgeshape : MyExplorer (wire, TopAbs_EDGE))
{
TopoDS_Edge edge = TopoDS::Edge(edgeshape);
curr++;
(*testout) << "edge nr " << curr << endl;
multithread.percent = 100 * curr / double (total);
if (multithread.terminate) return;
// TopoDS_Edge edge = TopoDS::Edge (exp3.Current());
if (BRep_Tool::Degenerated(edge))
{
//(*testout) << "ignoring degenerated edge" << endl;
continue;
}
if (!geom.emap.Contains(edge)) continue;
if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) ==
geom.vmap.FindIndex(TopExp::LastVertex (edge)))
{
GProp_GProps system;
BRepGProp::LinearProperties(edge, system);
if (system.Mass() < eps)
{
cout << "ignoring edge " << geom.emap.FindIndex (edge)
<< ". closed edge with length < " << eps << endl;
continue;
}
}
double s0, s1;
Handle(Geom2d_Curve) cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
int geomedgenr = geom.emap.FindIndex(edge);
Array<PointIndex> pnums;
Array<double> params;
// check for identifications
bool copy_identified = false;
if (auto it = geom.identifications.find(edge.TShape()); it != geom.identifications.end())
for (auto & ids : it->second)
{
cout << "edge has identification with trafo " << ids.name << ", inv = " << ids.inverse << endl;
int otherind = geom.emap.FindIndex(ids.other);
Array<Segment> othersegs;
for (auto & seg : mesh.LineSegments())
if (seg.edgenr == otherind)
othersegs.Append (seg);
if (othersegs.Size())
{
cout << "other has already segs" << endl;
copy_identified = true;
Array<PointIndex> pnums_other;
pnums_other.Append (othersegs[0][0]);
for (auto & seg : othersegs)
pnums_other.Append (seg[1]);
auto inv = ids.trafo.CalcInverse();
// for (auto & pi : pnums)
for (auto oi : Range(pnums_other))
{
PointIndex piother = pnums_other[pnums_other.Size()-oi-1];
Point<3> pother = mesh[piother];
Point<3> p = inv(pother);
bool found = false;
PointIndex pi;
for (PointIndex piv : vertexrange)
if (Dist2 (mesh[piv], p) < eps*eps)
{
pi = piv;
found = true;
}
if (!found)
pi = mesh.AddPoint (p);
// params.Add ( find parameter p );
double s0, s1;
Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, s0, s1);
GeomAPI_ProjectPointOnCurve proj(ng2occ(p), curve);
params.Append (proj.LowerDistanceParameter());
pnums.Append (pi);
mesh.GetIdentifications().Add (pi, piother, geomedgenr);
}
mesh.GetIdentifications().SetType(geomedgenr,Identifications::PERIODIC);
copy_identified = true;
break;
}
}
if (!copy_identified)
{
if (alledgepnums[geomedgenr-1].Size())
{
pnums = alledgepnums[geomedgenr-1];
params = alledgeparams[geomedgenr-1];
}
else
{
NgArray <MeshPoint> mp;
DivideEdge (edge, mp, params, mesh, mparam);
pnums.SetSize(mp.Size()+2);
if (!merge_solids)
{
pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;
pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;
}
else
{
Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));
Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));
pnums[0] = PointIndex::INVALID;
pnums.Last() = PointIndex::INVALID;
for (PointIndex pi : vertexrange)
{
if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;
if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;
}
}
for (size_t i = 1; i <= mp.Size(); i++)
{
bool exists = false;
tsearch.Start();
/*
// for (PointIndex j = first_ep; j < mesh.Points().End(); j++)
for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++)
if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
{
exists = true;
pnums[i] = j;
break;
}
*/
tsearch.Stop();
if (!exists)
pnums[i] = mesh.AddPoint (mp[i-1]);
}
}
alledgepnums[geomedgenr-1] = pnums;
alledgeparams[geomedgenr-1] = params;
}
auto name = geom.eprops[geom.emap.FindIndex(edge)-1]->name.value_or("");
mesh.SetCD2Name(geomedgenr, name);
(*testout) << "NP = " << mesh.GetNP() << endl;
//(*testout) << pnums[pnums.Size()-1] << endl;
double hpref = OCCGeometry::global_shape_properties[edge.TShape()].hpref;
// for (size_t i = 1; i <= mp.Size()+1; i++)
for (size_t i = 1; i < pnums.Size(); i++)
{
// edgenr++;
Segment seg;
seg[0] = pnums[i-1];
seg[1] = pnums[i];
// seg.edgenr = edgenr;
seg.edgenr = geomedgenr;
seg.si = facenr;
seg.epgeominfo[0].dist = params[i-1];
seg.epgeominfo[1].dist = params[i];
seg.epgeominfo[0].edgenr = geomedgenr;
seg.epgeominfo[1].edgenr = geomedgenr;
double s0 = params[i-1];
double s1 = params[i];
double delta = s1-s0;
s0 += 1e-10*delta; // fixes normal-vector roundoff problem when endpoint is cone-tip
s1 -= 1e-10*delta;
gp_Pnt2d p2d1 = cof->Value(s0);
gp_Pnt2d p2d2 = cof->Value(s1);
seg.epgeominfo[0].u = p2d1.X();
seg.epgeominfo[0].v = p2d1.Y();
seg.epgeominfo[1].u = p2d2.X();
seg.epgeominfo[1].v = p2d2.Y();
seg.singedge_left = hpref;
seg.singedge_right = hpref;
/*
if (occface->IsUPeriodic())
{
cout << "U Periodic" << endl;
if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) >
fabs(seg.epgeominfo[1].u-
(seg.epgeominfo[0].u-occface->UPeriod())))
seg.epgeominfo[0].u = p2d.X()+occface->UPeriod();
if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) >
fabs(seg.epgeominfo[1].u-
(seg.epgeominfo[0].u+occface->UPeriod())))
seg.epgeominfo[0].u = p2d.X()-occface->UPeriod();
}
if (occface->IsVPeriodic())
{
cout << "V Periodic" << endl;
if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) >
fabs(seg.epgeominfo[1].v-
(seg.epgeominfo[0].v-occface->VPeriod())))
seg.epgeominfo[0].v = p2d.Y()+occface->VPeriod();
if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) >
fabs(seg.epgeominfo[1].v-
(seg.epgeominfo[0].v+occface->VPeriod())))
seg.epgeominfo[0].v = p2d.Y()-occface->VPeriod();
}
*/
if (edge.Orientation() == TopAbs_REVERSED)
{
swap (seg[0], seg[1]);
swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
}
mesh.AddSegment (seg);
//edgesegments[geomedgenr-1]->Append(mesh.GetNSeg());
}
}
}
}
// for(i=1; i<=mesh.GetNSeg(); i++)
// (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si
// << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl;
// exit(10);
mesh.CalcSurfacesOfNode();
multithread.task = savetask;
}
void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh,
const MeshingParameters & mparam)
{
static Timer t("OCCMeshSurface"); RegionTimer r(t);
// int i, j, k;
// int changed;
const char * savetask = multithread.task;
multithread.task = "Surface meshing";
geom.facemeshstatus = 0;
int noldp = mesh.GetNP();
double starttime = GetTime();
NgArray<int, PointIndex::BASE> glob2loc(noldp);
//int projecttype = PARAMETERSPACE;
int projecttype = PARAMETERSPACE;
int notrys = 1;
int surfmesherror = 0;
for (int k = 1; k <= mesh.GetNFD(); k++)
bool OCCMeshFace (const OCCGeometry & geom, Mesh & mesh, FlatArray<int, PointIndex> glob2loc,
const MeshingParameters & mparam, int nr, int projecttype, bool delete_on_failure)
{
auto k = nr+1;
if(1==0 && !geom.fvispar[k-1].IsDrawable())
{
(*testout) << "ignoring face " << k << endl;
cout << "ignoring face " << k << endl;
continue;
return true;
}
// if(master_faces[k]!=k)
// continue;
(*testout) << "mesh face " << k << endl;
multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
geom.facemeshstatus[k-1] = -1;
FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
auto face = TopoDS::Face(geom.fmap(k));
auto fshape = face.TShape();
int oldnf = mesh.GetNSE();
@ -758,7 +262,7 @@ namespace netgen
static Timer tinit("init");
tinit.Start();
Meshing2OCCSurfaces meshing(geom, TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam);
Meshing2OCCSurfaces meshing(geom, face, bb, projecttype, mparam);
tinit.Stop();
@ -769,22 +273,22 @@ namespace netgen
else
PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)");
tprint.Stop();
if (surfmesherror)
cout << "Surface meshing error occurred before (in " << surfmesherror << " faces)" << endl;
// Meshing2OCCSurfaces meshing(f2, bb);
meshing.SetStartTime (starttime);
// meshing.SetStartTime (starttime);
//(*testout) << "Face " << k << endl << endl;
auto segments = geom.GetFace(k-1).GetBoundary(mesh);
if (meshing.GetProjectionType() == PLANESPACE)
{
static Timer t("MeshSurface: Find edges and points - Physical"); RegionTimer r(t);
int cntp = 0;
glob2loc = 0;
for (Segment & seg : mesh.LineSegments())
if (seg.si == k)
for (Segment & seg : segments)
// if (seg.si == k)
for (int j = 0; j < 2; j++)
{
PointIndex pi = seg[j];
@ -801,8 +305,9 @@ namespace netgen
{
Segment & seg = mesh.LineSegment(i);
*/
for (Segment & seg : mesh.LineSegments())
if (seg.si == k)
// for (Segment & seg : mesh.LineSegments())
for (Segment & seg : segments)
//if (seg.si == k)
{
PointGeomInfo gi0, gi1;
gi0.trignum = gi1.trignum = k;
@ -811,7 +316,9 @@ namespace netgen
gi1.u = seg.epgeominfo[1].u;
gi1.v = seg.epgeominfo[1].v;
//if(orientation & 1)
meshing.AddBoundaryElement (glob2loc[seg[0]], glob2loc[seg[1]], gi0, gi1);
}
}
else
@ -825,19 +332,21 @@ namespace netgen
if (mesh.LineSegment(i).si == k)
cntp+=2;
*/
for (Segment & seg : mesh.LineSegments())
if (seg.si == k)
cntp += 2;
cntp = 2*segments.Size();
//for (Segment & seg : mesh.LineSegments())
//if (seg.si == k)
//cntp += 2;
NgArray<PointGeomInfo> gis;
gis.SetAllocSize (cntp);
gis.SetSize (0);
for (int i = 1; i <= mesh.GetNSeg(); i++)
//for (int i = 1; i <= mesh.GetNSeg(); i++)
for(auto & seg : segments)
{
Segment & seg = mesh.LineSegment(i);
if (seg.si == k)
//Segment & seg = mesh.LineSegment(i);
//if (seg.si == k)
{
PointGeomInfo gi0, gi1;
gi0.trignum = gi1.trignum = k;
@ -898,8 +407,6 @@ namespace netgen
}
// Philippose - 15/01/2009
double maxh = min2(geom.face_maxh[k-1], OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].maxh);
//double maxh = mparam.maxh;
@ -939,168 +446,23 @@ namespace netgen
projecttype = PARAMETERSPACE;
static Timer t1("rest of loop"); RegionTimer reg1(t1);
if (res != MESHING2_OK)
{
if (notrys == 1)
bool meshing_failed = res != MESHING2_OK;
if(meshing_failed && delete_on_failure)
{
for (SurfaceElementIndex sei = noldsurfel; sei < mesh.GetNSE(); sei++)
mesh.Delete(sei);
mesh.Compress();
(*testout) << "retry Surface " << k << endl;
k--;
// projecttype*=-1;
projecttype = PLANESPACE;
notrys++;
continue;
}
else
{
geom.facemeshstatus[k-1] = -1;
PrintError ("Problem in Surface mesh generation");
surfmesherror++;
// throw NgException ("Problem in Surface mesh generation");
}
}
else
{
geom.facemeshstatus[k-1] = 1;
}
notrys = 1;
for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++)
mesh[sei].SetIndex (k);
auto n_illegal_trigs = mesh.FindIllegalTrigs();
PrintMessage (3, n_illegal_trigs, " illegal triangles");
return meshing_failed;
}
// ofstream problemfile("occmesh.rep");
// problemfile << "SURFACEMESHING" << endl << endl;
if (surfmesherror)
{
cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl;
cout << "SURFACE MESHING ERROR OCCURRED IN " << surfmesherror << " FACES:" << endl;
for (int i = 1; i <= geom.fmap.Extent(); i++)
if (geom.facemeshstatus[i-1] == -1)
{
cout << "Face " << i << endl;
// problemfile << "problem with face " << i << endl;
// problemfile << "vertices: " << endl;
TopExp_Explorer exp0,exp1,exp2;
for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() )
{
TopoDS_Wire wire = TopoDS::Wire(exp0.Current());
for ( exp1.Init(wire,TopAbs_EDGE); exp1.More(); exp1.Next() )
{
TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
for ( exp2.Init(edge,TopAbs_VERTEX); exp2.More(); exp2.Next() )
{
TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current());
gp_Pnt point = BRep_Tool::Pnt(vertex);
// problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl;
}
}
}
// problemfile << endl;
}
cout << endl << endl;
cout << "for more information open IGES/STEP Topology Explorer" << endl;
// problemfile.close();
throw NgException ("Problem in Surface mesh generation");
}
else
{
// problemfile << "OK" << endl << endl;
// problemfile.close();
}
for (int i = 0; i < mesh.GetNFD(); i++)
mesh.SetBCName (i, mesh.GetFaceDescriptor(i+1).GetBCName());
multithread.task = savetask;
}
void OCCOptimizeSurface(OCCGeometry & geom, Mesh & mesh,
const MeshingParameters & mparam)
{
const char * savetask = multithread.task;
multithread.task = "Optimizing surface";
static Timer timer_opt2d("Optimization 2D");
timer_opt2d.Start();
for (int k = 1; k <= mesh.GetNFD(); k++)
{
// if (k != 42) continue;
// if (k != 36) continue;
// (*testout) << "optimize face " << k << endl;
multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
PrintMessage (1, "Optimize Surface ", k);
for (int i = 1; i <= mparam.optsteps2d; i++)
{
// (*testout) << "optstep " << i << endl;
if (multithread.terminate) return;
{
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.EdgeSwapping (i > mparam.optsteps2d/2);
}
if (multithread.terminate) return;
{
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mparam);
}
{
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.CombineImprove ();
}
if (multithread.terminate) return;
{
MeshOptimize2d meshopt(mesh);
meshopt.SetFaceIndex (k);
meshopt.SetImproveEdges (0);
meshopt.SetMetricWeight (mparam.elsizeweight);
meshopt.SetWriteStatus (0);
meshopt.ImproveMesh (mparam);
}
}
}
mesh.CalcSurfacesOfNode();
mesh.Compress();
timer_opt2d.Stop();
multithread.task = savetask;
}
void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh,
const MeshingParameters & mparam, const OCCParameters& occparam)

688
libsrc/occ/occgeom.cpp
View File

@ -1,52 +1,69 @@
#ifdef OCCGEOMETRY
#include <mystdlib.h>
#include <occgeom.hpp>
#include <core/register_archive.hpp>
#include <cstdio>
#include "ShapeAnalysis_ShapeTolerance.hxx"
#include "ShapeAnalysis_ShapeContents.hxx"
#include "ShapeAnalysis_CheckSmallFace.hxx"
#include "ShapeAnalysis_DataMapOfShapeListOfReal.hxx"
#include "ShapeAnalysis_Surface.hxx"
#include <set>
#include "BRepCheck_Analyzer.hxx"
#include "BRepLib.hxx"
#include "ShapeBuild_ReShape.hxx"
#include "ShapeFix.hxx"
#include "ShapeFix_FixSmallFace.hxx"
#include "Partition_Spliter.hxx"
#include "BRepAlgoAPI_Fuse.hxx"
#include "Interface_InterfaceModel.hxx"
#include <mystdlib.h>
#include <core/register_archive.hpp>
#include "XSControl_WorkSession.hxx"
#include "XSControl_TransferReader.hxx"
#include "StepRepr_RepresentationItem.hxx"
#include "StepBasic_ProductDefinitionRelationship.hxx"
#include "Transfer_TransientProcess.hxx"
#include "TransferBRep.hxx"
#ifndef _Standard_Version_HeaderFile
#include <Standard_Version.hxx>
#endif
#include "occ_vertex.hpp"
#include "occ_edge.hpp"
#include "occ_face.hpp"
#include "occ_solid.hpp"
#include "occgeom.hpp"
#include <BOPAlgo_Builder.hxx>
#include <BRepBndLib.hxx>
#include <BRepBuilderAPI_MakeSolid.hxx>
#include <BRepBuilderAPI_MakeVertex.hxx>
#include <BRepCheck_Analyzer.hxx>
#include <BRepExtrema_DistShapeShape.hxx>
#include <BRepGProp.hxx>
#include <BRepLib.hxx>
#include <BRepMesh_IncrementalMesh.hxx>
#include <BRepOffsetAPI_Sewing.hxx>
#include <BRepTools.hxx>
#include <IGESCAFControl_Reader.hxx>
#include <IGESControl_Writer.hxx>
#include <Interface_InterfaceModel.hxx>
#include <Interface_Static.hxx>
#include <Partition_Spliter.hxx>
#include <STEPCAFControl_Writer.hxx>
#include <STEPConstruct.hxx>
#include <Transfer_FinderProcess.hxx>
#include <TDataStd_Name.hxx>
#include <XCAFPrs.hxx>
#include <STEPControl_Writer.hxx>
#include <ShapeAnalysis_CheckSmallFace.hxx>
#include <ShapeAnalysis_DataMapOfShapeListOfReal.hxx>
#include <ShapeAnalysis_ShapeContents.hxx>
#include <ShapeBuild_ReShape.hxx>
#include <ShapeFix_Face.hxx>
#include <ShapeFix_FixSmallFace.hxx>
#include <ShapeFix_Shape.hxx>
#include <ShapeFix_Wire.hxx>
#include <ShapeFix_Wireframe.hxx>
#include <StepBasic_ProductDefinitionRelationship.hxx>
#include <StepRepr_RepresentationItem.hxx>
#include <StlAPI_Writer.hxx>
#include <TopoDS_Shape.hxx>
#include <TransferBRep.hxx>
#include <Transfer_FinderProcess.hxx>
#include <Transfer_TransientProcess.hxx>
#include <XCAFApp_Application.hxx>
#include <XCAFDoc_ColorTool.hxx>
#include <XCAFDoc_DocumentTool.hxx>
#include <XCAFDoc_ShapeTool.hxx>
#include <XCAFPrs.hxx>
#include <XCAFPrs_Style.hxx>
#include <TopTools_ShapeMapHasher.hxx>
#include <NCollection_IndexedDataMap.hxx>
#include <StepShape_TopologicalRepresentationItem.hxx>
#include <XSControl_TransferReader.hxx>
#include <XSControl_WorkSession.hxx>
#if OCC_VERSION_HEX < 0x070000
// pass
#elif OCC_VERSION_HEX < 0x070200
#include "StlTransfer.hxx"
#include "TopoDS_Iterator.hxx"
#include <StlTransfer.hxx>
#include <TopoDS_Iterator.hxx>
#else
#include "TopoDS_Iterator.hxx"
#include <TopoDS_Iterator.hxx>
#endif
namespace netgen
@ -57,6 +74,67 @@ namespace netgen
std::map<Handle(TopoDS_TShape), ShapeProperties> OCCGeometry::global_shape_properties;
std::map<Handle(TopoDS_TShape), std::vector<OCCIdentification>> OCCGeometry::identifications;
TopoDS_Shape ListOfShapes::Max(gp_Vec dir)
{
double maxval = -1e99;
TopoDS_Shape maxshape;
for (auto shape : *this)
{
GProp_GProps props;
gp_Pnt center;
switch (shape.ShapeType())
{
case TopAbs_VERTEX:
center = BRep_Tool::Pnt (TopoDS::Vertex(shape)); break;
case TopAbs_FACE:
BRepGProp::SurfaceProperties (shape, props);
center = props.CentreOfMass();
break;
default:
BRepGProp::LinearProperties(shape, props);
center = props.CentreOfMass();
}
double val = center.X()*dir.X() + center.Y()*dir.Y() + center.Z() * dir.Z();
if (val > maxval)
{
maxval = val;
maxshape = shape;
}
}
return maxshape;
}
TopoDS_Shape ListOfShapes::Nearest(gp_Pnt pnt)
{
double mindist = 1e99;
TopoDS_Shape nearestshape;
auto vertex = BRepBuilderAPI_MakeVertex (pnt).Vertex();
for (auto shape : *this)
{
double dist = BRepExtrema_DistShapeShape(shape, vertex).Value();
if (dist < mindist)
{
nearestshape = shape;
mindist = dist;
}
}
return nearestshape;
}
ListOfShapes ListOfShapes::SubShapes(TopAbs_ShapeEnum type) const
{
std::set<TopoDS_Shape, ShapeLess> unique_shapes;
for(const auto& shape : *this)
for(TopExp_Explorer e(shape, type); e.More(); e.Next())
unique_shapes.insert(e.Current());
ListOfShapes sub;
for(const auto& shape : unique_shapes)
sub.push_back(shape);
return sub;
}
OCCGeometry::OCCGeometry(const TopoDS_Shape& _shape, int aoccdim, bool copy)
{
if(copy)
@ -64,14 +142,14 @@ namespace netgen
auto filename = GetTempFilename();
step_utils::WriteSTEP(_shape, filename.c_str());
LoadOCCInto(this, filename.c_str());
occdim = aoccdim;
dimension = aoccdim;
std::remove(filename.c_str());
}
else
{
shape = _shape;
changed = 1;
occdim = aoccdim;
dimension = aoccdim;
BuildFMap();
CalcBoundingBox();
PrintContents (this);
@ -117,23 +195,23 @@ namespace netgen
OCCSetLocalMeshSize(*this, mesh, mparam, occparam);
}
void OCCGeometry :: FindEdges(Mesh& mesh,
const MeshingParameters& mparam) const
bool OCCGeometry :: MeshFace(Mesh& mesh,
const MeshingParameters& mparam, int nr, FlatArray<int, PointIndex> glob2loc) const
{
OCCFindEdges(*this, mesh, mparam);
}
bool failed = OCCMeshFace(*this, mesh, glob2loc, mparam, nr, PLANESPACE, true);
if(failed)
failed = OCCMeshFace(*this, mesh, glob2loc, mparam, nr, PARAMETERSPACE, false);
void OCCGeometry :: MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) const
if(failed)
{
OCCMeshSurface(*this, mesh, mparam);
facemeshstatus[nr] = -1;
PrintError ("Problem in Surface mesh generation");
}
void OCCGeometry :: FinalizeMesh(Mesh& mesh) const
else
{
for (int i = 0; i < mesh.GetNDomains(); i++)
if (auto name = sprops[i]->name)
mesh.SetMaterial (i+1, *name);
facemeshstatus[nr] = 1;
}
return failed;
}
void OCCGeometry :: PrintNrShapes ()
@ -1029,31 +1107,97 @@ namespace netgen
fsingular = esingular = vsingular = false;
sprops.SetSize(somap.Extent());
for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
// Add shapes
for(auto v : GetVertices(shape))
{
auto s = e.Current();
sprops[somap.FindIndex(s)-1] = &global_shape_properties[s.TShape()];
auto tshape = v.TShape();
if(vertex_map.count(tshape)!=0)
continue;
vertex_map[tshape] = vertices.Size();
auto occ_vertex = make_unique<OCCVertex>(TopoDS::Vertex(v));
if(global_shape_properties.count(tshape)>0)
occ_vertex->properties = global_shape_properties[tshape];
vertices.Append(std::move(occ_vertex));
}
fprops.SetSize(fmap.Extent());
for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
for(auto e : GetEdges(shape))
{
auto s = e.Current();
fprops[fmap.FindIndex(s)-1] = &global_shape_properties[s.TShape()];
auto tshape = e.TShape();
auto edge = TopoDS::Edge(e);
if(edge_map.count(tshape)!=0)
continue;
if(BRep_Tool::Degenerated(edge))
continue;
edge_map[tshape] = edges.Size();
auto occ_edge = make_unique<OCCEdge>(edge);
occ_edge->properties = global_shape_properties[tshape];
edges.Append(std::move(occ_edge));
}
eprops.SetSize(emap.Extent());
/*
for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next())
for(auto f : GetFaces(shape))
{
auto s = e.Current();
eprops[emap.FindIndex(s)-1] = &global_shape_properties[s.TShape()];
auto tshape = f.TShape();
if(face_map.count(tshape)==0)
{
auto k = faces.Size();
face_map[tshape] = k;
auto occ_face = make_unique<OCCFace>(f);
if(global_shape_properties.count(tshape)>0)
occ_face->properties = global_shape_properties[tshape];
faces.Append(std::move(occ_face));
}
*/
for (auto [nr,s] : Enumerate(emap))
eprops[nr-1] = &global_shape_properties[s.TShape()];
}
for(auto s : GetSolids(shape))
{
auto tshape = s.TShape();
int k;
if(solid_map.count(tshape)==0)
{
k = solids.Size();
solid_map[tshape] = k;
auto occ_solid = make_unique<OCCSolid>(s);
if(global_shape_properties.count(tshape)>0)
occ_solid->properties = global_shape_properties[tshape];
solids.Append(std::move(occ_solid));
}
for(auto f : GetFaces(s))
{
auto face_nr = face_map[f.TShape()];
auto & face = faces[face_nr];
if(face->domin==-1)
face->domin = k;
else
face->domout = k;
}
}
// Add identifications
auto add_identifications = [&](auto & shapes, auto & shape_map)
{
for(auto &[tshape, nr] : shape_map)
if(identifications.count(tshape))
for(auto & ident : identifications[tshape])
{
ShapeIdentification si{
shapes[shape_map[ident.from]].get(),
shapes[shape_map[ident.to]].get(),
ident.trafo,
ident.type,
ident.name
};
shapes[nr]->identifications.Append(si);
}
};
add_identifications( vertices, vertex_map );
add_identifications( edges, edge_map );
add_identifications( faces, face_map );
ProcessIdentifications();
bounding_box = ::netgen::GetBoundingBox( shape );
}
@ -1156,286 +1300,11 @@ namespace netgen
void OCCGeometry :: CalcBoundingBox ()
{
Bnd_Box bb;
#if OCC_VERSION_HEX < 0x070000
BRepBndLib::Add (shape, bb);
#else
BRepBndLib::Add ((const TopoDS_Shape) shape, bb,(Standard_Boolean)true);
#endif
double x1,y1,z1,x2,y2,z2;
bb.Get (x1,y1,z1,x2,y2,z2);
Point<3> p1 = Point<3> (x1,y1,z1);
Point<3> p2 = Point<3> (x2,y2,z2);
(*testout) << "Bounding Box = [" << p1 << " - " << p2 << "]" << endl;
boundingbox = Box<3> (p1,p2);
boundingbox = ::netgen::GetBoundingBox(shape);
(*testout) << "Bounding Box = [" << boundingbox.PMin() << " - " << boundingbox.PMax() << "]" << endl;
SetCenter();
}
PointGeomInfo OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const
{
static int cnt = 0;
if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
gp_Pnt pnt(p(0), p(1), p(2));
double u,v;
Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) );
suval.Coord( u, v);
pnt = thesurf->Value( u, v );
PointGeomInfo gi;
gi.trignum = surfi;
gi.u = u;
gi.v = v;
p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
return gi;
}
bool OCCGeometry :: ProjectPointGI(int surfind, Point<3>& p, PointGeomInfo& gi) const
{
double u = gi.u;
double v = gi.v;
Point<3> hp = p;
if (FastProject (surfind, hp, u, v))
{
p = hp;
return 1;
}
ProjectPoint (surfind, p);
return CalcPointGeomInfo (surfind, gi, p);
}
void OCCGeometry :: ProjectPointEdge(int surfind, INDEX surfind2,
Point<3> & p, EdgePointGeomInfo* gi) const
{
TopExp_Explorer exp0, exp1;
bool done = false;
Handle(Geom_Curve) c;
for (exp0.Init(fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next())
for (exp1.Init(fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next())
{
if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current())))
{
done = true;
double s0, s1;
c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1);
}
}
gp_Pnt pnt(p(0), p(1), p(2));
GeomAPI_ProjectPointOnCurve proj(pnt, c);
pnt = proj.NearestPoint();
p(0) = pnt.X();
p(1) = pnt.Y();
p(2) = pnt.Z();
}
bool OCCGeometry :: FastProject (int surfi, Point<3> & ap, double& u, double& v) const
{
gp_Pnt p(ap(0), ap(1), ap(2));
Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
gp_Pnt x = surface->Value (u,v);
if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
gp_Vec du, dv;
surface->D1(u,v,x,du,dv);
int count = 0;
gp_Pnt xold;
gp_Vec n;
double det, lambda, mu;
do {
count++;
n = du^dv;
det = Det3 (n.X(), du.X(), dv.X(),
n.Y(), du.Y(), dv.Y(),
n.Z(), du.Z(), dv.Z());
if (det < 1e-15) return false;
lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
n.Y(), p.Y()-x.Y(), dv.Y(),
n.Z(), p.Z()-x.Z(), dv.Z())/det;
mu = Det3 (n.X(), du.X(), p.X()-x.X(),
n.Y(), du.Y(), p.Y()-x.Y(),
n.Z(), du.Z(), p.Z()-x.Z())/det;
u += lambda;
v += mu;
xold = x;
surface->D1(u,v,x,du,dv);
} while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
// (*testout) << "FastProject count: " << count << endl;
if (count == 50) return false;
ap = Point<3> (x.X(), x.Y(), x.Z());
return true;
}
Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* geominfo) const
{
if(geominfo)
{
gp_Pnt pnt;
gp_Vec du, dv;
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
occface->D1(geominfo->u,geominfo->v,pnt,du,dv);
auto n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
Vec<3>(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
return n;
}
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) );
suval.Coord( u, v);
pnt = occface->Value( u, v );
gp_Vec du, dv;
occface->D1(u,v,pnt,du,dv);
/*
if (!occface->IsCNu (1) || !occface->IsCNv (1))
(*testout) << "SurfOpt: Differentiation FAIL" << endl;
*/
auto n = Cross (Vec3d(du.X(), du.Y(), du.Z()),
Vec3d(dv.X(), dv.Y(), dv.Z()));
n.Normalize();
if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
return n;
}
bool OCCGeometry :: CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const
{
Standard_Real u,v;
gp_Pnt pnt(p(0), p(1), p(2));
Handle(Geom_Surface) occface;
occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
/*
GeomAPI_ProjectPointOnSurf proj(pnt, occface);
if (proj.NbPoints() < 1)
{
cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
<< endl;
cout << p << endl;
return 0;
}
proj.LowerDistanceParameters (u, v);
*/
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) );
suval.Coord( u, v);
//pnt = occface->Value( u, v );
gi.u = u;
gi.v = v;
return true;
}
void OCCGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const
{
Point<3> hnewp;
hnewp = p1+secpoint*(p2-p1);
if (surfi > 0)
{
double u = gi1.u+secpoint*(gi2.u-gi1.u);
double v = gi1.v+secpoint*(gi2.v-gi1.v);
auto savept = hnewp;
if (!FastProject(surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2))
{
// cout << "Fast projection to surface fails! Using OCC projection" << endl;
hnewp = savept;
ProjectPoint(surfi, hnewp);
}
newgi.trignum = 1;
newgi.u = u;
newgi.v = v;
}
newp = hnewp;
}
void OCCGeometry :: PointBetweenEdge(const Point<3> & p1,
const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const
{
double s0, s1;
newp = p1+secpoint*(p2-p1);
if(ap1.edgenr > emap.Size() || ap1.edgenr == 0)
return;
gp_Pnt pnt(newp(0), newp(1), newp(2));
GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(emap(ap1.edgenr)), s0, s1));
pnt = proj.NearestPoint();
newp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
newgi = ap1;
};
// void OCCGeometry :: WriteOCC_STL(char * filename)
// {
@ -1699,26 +1568,18 @@ namespace netgen
auto occ_hash = key.HashCode(1<<31UL);
return std::hash<decltype(occ_hash)>()(occ_hash);
};
std::unordered_map<TopoDS_Shape, int, decltype(my_hash)> shape_map(10, my_hash);
Array<TopoDS_Shape> shape_list;
std::map<Handle(TopoDS_TShape), int> tshape_map;
Array<Handle(TopoDS_TShape)> tshape_list;
ar & occdim;
ar & dimension;
for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE })
for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
{
auto ds = e.Current();
auto ts = ds.TShape();
if(shape_map.count(ds)==0)
{
shape_map[ds] = shape_list.Size();
shape_list.Append(ds);
}
if(tshape_map.count(ts)==0)
{
tshape_map[ts] = shape_list.Size();
tshape_map[ts] = tshape_list.Size();
tshape_list.Append(ts);
}
}
@ -1741,12 +1602,18 @@ namespace netgen
for(auto i : Range(n_idents))
{
auto & id = idents[i];
int shape_id;
int id_from, id_to;
if(ar.Output())
shape_id = shape_map[id.other];
ar & shape_id & id.trafo & id.inverse & id.name;
{
id_from = tshape_map[id.from];
id_to = tshape_map[id.to];
}
ar & id_from & id_to & id.trafo & id.name;
if(ar.Input())
id.other = shape_list[shape_id];
{
id.from = tshape_list[id_from];
id.to = tshape_list[id_to];
}
}
}
}
@ -2053,6 +1920,97 @@ namespace netgen
return false;
}
Point<3> GetCenter(const TopoDS_Shape & shape)
{
GProp_GProps props;
BRepGProp::LinearProperties(shape, props);
return occ2ng( props.CentreOfMass() );
}
void OCCGeometry :: IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type)
{
auto cme = GetCenter(me);
auto cyou = GetCenter(you);
Transformation<3> trafo{cyou-cme};
identifications[me.TShape()].push_back( {me.TShape(), you.TShape(), Transformation<3>(cyou - cme), name, type} );
auto vme = GetVertices(me);
auto vyou = GetVertices(you);
Point<3> pme0 = trafo(occ2ng(vme[0]));
Point<3> pme1 = trafo(occ2ng(vme[1]));
Point<3> pyou = occ2ng(vyou[0]);
bool do_swap = Dist(pme1, pyou) < Dist(pme0, pyou);
for(auto i : Range(2))
identifications[vme[i].TShape()].push_back( {vme[i].TShape(), vyou[do_swap ? 1-i : i].TShape(), trafo, name, type} );
}
bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you)
{
if(me.ShapeType() != you.ShapeType()) return false;
Bnd_Box bbox;
BRepBndLib::Add(me, bbox);
BRepBndLib::Add(you, bbox);
BoxTree<3> tree( occ2ng(bbox.CornerMin()), occ2ng(bbox.CornerMax()) );
Point<3> c_me = GetCenter(me);
Point<3> c_you = GetCenter(you);
if(tree.GetTolerance() < Dist(trafo(c_me), c_you))
return false;
std::map<T_Shape, T_Shape> vmap;
auto verts_me = GetVertices(me);
for (auto i : Range(verts_me.size()))
{
auto s = verts_me[i].TShape();
if(vmap.count(s)>0)
throw Exception("vertex mapped twice!");
auto p = trafo(occ2ng(s));
tree.Insert( p, i );
vmap[s] = nullptr;
}
bool all_verts_mapped = true;
for (auto vert : GetVertices(you))
{
auto s = vert.TShape();
auto p = occ2ng(s);
bool vert_mapped = false;
tree.GetFirstIntersecting( p, p, [&](size_t i ) {
vmap[verts_me[i].TShape()] = s;
vert_mapped = true;
return true;
});
if(!vert_mapped)
{
all_verts_mapped = false;
break;
}
}
return all_verts_mapped;
}
void OCCGeometry :: IdentifyFaces(const TopoDS_Shape & solid, const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type)
{
auto cme = GetCenter(me);
auto cyou = GetCenter(you);
Transformation<3> trafo(cyou-cme);
identifications[me.TShape()].push_back
(OCCIdentification { me.TShape(), you.TShape(), trafo, name, type });
auto edges_me = GetEdges(me);
auto edges_you = GetEdges(you);
for (auto e_me : edges_me)
for (auto e_you : edges_you)
if(IsMappedShape(trafo, e_me, e_you))
IdentifyEdges(e_me, e_you, name, type);
}
void OCCParameters :: Print(ostream & ost) const
{
ost << "OCC Parameters:" << endl
@ -2063,6 +2021,18 @@ namespace netgen
DLL_HEADER extern OCCParameters occparam;
OCCParameters occparam;
// int OCCGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
// {
// return OCCGenerateMesh (*this, mesh, mparam, occparam);
@ -2208,8 +2178,7 @@ namespace netgen
for(auto & ident : identifications)
{
Array<Handle(StepRepr_RepresentationItem)> items;
items.Append(STEPConstruct::FindEntity(finder, ident.other));
items.Append(MakeInt(static_cast<int>(ident.inverse)));
// items.Append(STEPConstruct::FindEntity(finder, ident.other)); // TODO!
auto & m = ident.trafo.GetMatrix();
for(auto i : Range(9))
items.Append(MakeReal(m(i)));
@ -2239,8 +2208,7 @@ namespace netgen
auto id_item = Handle(StepRepr_CompoundRepresentationItem)::DownCast(idents->ItemElementValue(i));
OCCIdentification ident;
ident.name = id_item->Name()->ToCString();
ident.other = TransferBRep::ShapeResult(transProc->Find(id_item->ItemElementValue(1)));
ident.inverse = static_cast<bool>(ReadInt(id_item->ItemElementValue(2)));
// ident.other = TransferBRep::ShapeResult(transProc->Find(id_item->ItemElementValue(1))); /TODO!
auto & m = ident.trafo.GetMatrix();
for(auto i : Range(9))

254
libsrc/occ/occgeom.hpp
View File

@ -10,99 +10,27 @@
#ifdef OCCGEOMETRY
#include <meshing.hpp>
#include "occ_utils.hpp"
#include "occmeshsurf.hpp"
#include <Standard_Version.hxx>
#include "BRep_Tool.hxx"
#include "Geom_Curve.hxx"
#include "Geom2d_Curve.hxx"
#include "Geom_Surface.hxx"
#include "GeomAPI_ProjectPointOnSurf.hxx"
#include "GeomAPI_ProjectPointOnCurve.hxx"
#include "BRepTools.hxx"
#include "TopExp.hxx"
#include "BRepBuilderAPI_MakeVertex.hxx"
#include "BRepBuilderAPI_MakeShell.hxx"
#include "BRepBuilderAPI_MakeSolid.hxx"
#include "BRepOffsetAPI_Sewing.hxx"
#include "BRepLProp_SLProps.hxx"
#include "BRepAdaptor_Surface.hxx"
#include "Poly_Triangulation.hxx"
#include "Poly_Array1OfTriangle.hxx"
#include "TColgp_Array1OfPnt2d.hxx"
#include "Poly_Triangle.hxx"
#include "GProp_GProps.hxx"
#include "BRepGProp.hxx"
#include "gp_Pnt.hxx"
#include "TopoDS.hxx"
#include "TopoDS_Solid.hxx"
#include "TopExp_Explorer.hxx"
#include "TopTools_ListIteratorOfListOfShape.hxx"
#include "TopoDS_Wire.hxx"
#include "BRepTools_WireExplorer.hxx"
#include "TopTools_IndexedMapOfShape.hxx"
#include "BRepLProp_CLProps.hxx"
#include "BRepAdaptor_Curve.hxx"
#include "TopoDS_Shape.hxx"
#include "TopoDS_Face.hxx"
#include "IGESToBRep_Reader.hxx"
#include "Interface_Static.hxx"
#include "GeomAPI_ExtremaCurveCurve.hxx"
#include "Standard_ErrorHandler.hxx"
#include "Standard_Failure.hxx"
#include "ShapeUpgrade_ShellSewing.hxx"
#include "ShapeFix_Shape.hxx"
#include "ShapeFix_Wireframe.hxx"
#include "BRepMesh_IncrementalMesh.hxx"
#include "BRepBndLib.hxx"
#include "Bnd_Box.hxx"
#include "ShapeAnalysis.hxx"
#include "ShapeBuild_ReShape.hxx"
#include "BOPAlgo_Builder.hxx"
// Philippose - 29/01/2009
// OpenCascade XDE Support
// Include support for OpenCascade XDE Features
#include "TDocStd_Document.hxx"
#include "Quantity_Color.hxx"
#include "XCAFApp_Application.hxx"
#include "XCAFDoc_ShapeTool.hxx"
#include "XCAFDoc_Color.hxx"
#include "XCAFDoc_ColorTool.hxx"
#include "XCAFDoc_ColorType.hxx"
#include "XCAFDoc_LayerTool.hxx"
#include "XCAFDoc_DimTolTool.hxx"
#include "XCAFDoc_MaterialTool.hxx"
#include "XCAFDoc_DocumentTool.hxx"
#include "TDF_Label.hxx"
#include "TDF_LabelSequence.hxx"
#include "STEPCAFControl_Reader.hxx"
#include "STEPCAFControl_Writer.hxx"
#include "IGESCAFControl_Reader.hxx"
#include "IGESCAFControl_Writer.hxx"
#include "IGESControl_Reader.hxx"
#include "STEPControl_Reader.hxx"
#include "IGESControl_Writer.hxx"
#include "STEPControl_Writer.hxx"
#include <StepRepr_ValueRepresentationItem.hxx>
#include <StepRepr_IntegerRepresentationItem.hxx>
#include <StepRepr_CompoundRepresentationItem.hxx>
#include <Quantity_ColorRGBA.hxx>
#include <STEPCAFControl_Reader.hxx>
#include <StepBasic_MeasureValueMember.hxx>
#include "StlAPI_Writer.hxx"
#include "STEPControl_StepModelType.hxx"
#include <StepRepr_CompoundRepresentationItem.hxx>
#include <StepRepr_IntegerRepresentationItem.hxx>
#include <StepRepr_ValueRepresentationItem.hxx>
#include <TCollection_HAsciiString.hxx>
#include <TDocStd_Document.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Shape.hxx>
#include <Transfer_FinderProcess.hxx>
#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
#define OCC_HAVE_HISTORY
#endif
namespace netgen
{
#include "occmeshsurf.hpp"
// extern DLL_HEADER MeshingParameters mparam;
@ -117,28 +45,6 @@ namespace netgen
#define OCCGEOMETRYVISUALIZATIONHALFCHANGE 2 // Redraw
inline Point<3> occ2ng (const gp_Pnt & p)
{
return Point<3> (p.X(), p.Y(), p.Z());
}
inline Point<2> occ2ng (const gp_Pnt2d & p)
{
return Point<2> (p.X(), p.Y());
}
inline Vec<3> occ2ng (const gp_Vec & v)
{
return Vec<3> (v.X(), v.Y(), v.Z());
}
inline gp_Pnt ng2occ (const Point<3> & p)
{
return gp_Pnt(p(0), p(1), p(2));
}
class EntityVisualizationCode
{
int code;
@ -219,108 +125,20 @@ namespace netgen
};
class ShapeProperties
{
public:
optional<string> name;
optional<Vec<4>> col;
double maxh = 1e99;
double hpref = 0; // number of hp refinement levels (will be multiplied by factor later)
void Merge(const ShapeProperties & prop2)
{
if (prop2.name) name = prop2.name;
if (prop2.col) col = prop2.col;
maxh = min2(maxh, prop2.maxh);
}
void DoArchive(Archive& ar)
{
ar & name & col & maxh & hpref;
}
};
class OCCIdentification
{
public:
TopoDS_Shape other;
Transformation<3> trafo;
bool inverse;
string name;
};
class MyExplorer
{
class Iterator
{
TopExp_Explorer exp;
public:
Iterator (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid)
: exp(ashape, atoFind, atoAvoid) { }
auto operator*() { return exp.Current(); }
Iterator & operator++() { exp.Next(); return *this; }
bool operator!= (nullptr_t nu) { return exp.More(); }
};
public:
TopoDS_Shape shape;
TopAbs_ShapeEnum toFind;
TopAbs_ShapeEnum toAvoid;
MyExplorer (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid = TopAbs_SHAPE)
: shape(ashape), toFind(atoFind), toAvoid(atoAvoid) { ; }
Iterator begin() { return Iterator(shape, toFind, toAvoid); }
auto end() { return nullptr; }
};
inline auto Explore (TopoDS_Shape shape, TopAbs_ShapeEnum toFind, TopAbs_ShapeEnum toAvoid = TopAbs_SHAPE)
{
return MyExplorer (shape, toFind, toAvoid);
}
class IndexMapIterator
{
class Iterator
{
const TopTools_IndexedMapOfShape & indmap;
int i;
public:
Iterator (const TopTools_IndexedMapOfShape & aindmap, int ai)
: indmap(aindmap), i(ai) { ; }
auto operator*() { return tuple(i, indmap(i)); }
Iterator & operator++() { i++; return *this; }
bool operator!= (const Iterator & i2) { return i != i2.i; }
};
public:
const TopTools_IndexedMapOfShape & indmap;
IndexMapIterator (const TopTools_IndexedMapOfShape & aindmap) : indmap(aindmap) { }
Iterator begin() { return Iterator(indmap, 1); }
Iterator end() { return Iterator(indmap, indmap.Extent()+1); }
};
inline auto Enumerate (const TopTools_IndexedMapOfShape & indmap)
{
return IndexMapIterator(indmap);
}
class DLL_HEADER OCCGeometry : public NetgenGeometry
{
Point<3> center;
OCCParameters occparam;
public:
static std::map<Handle(TopoDS_TShape), ShapeProperties> global_shape_properties;
static std::map<Handle(TopoDS_TShape), std::vector<OCCIdentification>> identifications;
static std::map<T_Shape, ShapeProperties> global_shape_properties;
static std::map<T_Shape, std::vector<OCCIdentification>> identifications;
TopoDS_Shape shape;
TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap;
TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; // legacy maps
NgArray<bool> fsingular, esingular, vsingular;
Box<3> boundingbox;
// should we use 1-based arrays (JS->MH) ?
Array<ShapeProperties*> fprops, eprops, sprops; // pointers to the gobal property map
std::map<T_Shape, int> edge_map, vertex_map, face_map, solid_map;
mutable int changed;
mutable NgArray<int> facemeshstatus;
@ -350,8 +168,6 @@ namespace netgen
bool makesolids;
bool splitpartitions;
int occdim = 3; // meshing is always done 3D, changed to 2D later of occdim=2
OCCGeometry()
{
somap.Clear();
@ -372,36 +188,15 @@ namespace netgen
void Analyse(Mesh& mesh,
const MeshingParameters& mparam) const override;
void FindEdges(Mesh& mesh,
const MeshingParameters& mparam) const override;
void MeshSurface(Mesh& mesh,
const MeshingParameters& mparam) const override;
void FinalizeMesh(Mesh& mesh) const override;
bool MeshFace(Mesh& mesh, const MeshingParameters& mparam,
int nr, FlatArray<int, PointIndex> glob2loc) const override;
// void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const override {}
void Save (string filename) const override;
void SaveToMeshFile (ostream & /* ost */) const override;
void DoArchive(Archive& ar) override;
PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override;
void ProjectPointEdge (int surfind, int surfind2, Point<3> & p,
EdgePointGeomInfo* gi = nullptr) const override;
bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override;
bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi1, int surfi2,
const EdgePointGeomInfo & ap1,
const EdgePointGeomInfo & ap2,
Point<3> & newp, EdgePointGeomInfo & newgi) const override;
void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
int surfi,
const PointGeomInfo & gi1,
const PointGeomInfo & gi2,
Point<3> & newp, PointGeomInfo & newgi) const override;
void BuildFMap();
auto GetShape() const { return shape; }
@ -546,9 +341,11 @@ namespace netgen
bool ErrorInSurfaceMeshing ();
// void WriteOCC_STL(char * filename);
static void IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type);
static void IdentifyFaces(const TopoDS_Shape & solid,const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type);
private:
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
//bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
};
@ -564,11 +361,8 @@ namespace netgen
DLL_HEADER extern void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam,
const OCCParameters& occparam);
DLL_HEADER extern void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
DLL_HEADER extern void OCCOptimizeSurface (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
DLL_HEADER extern void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
DLL_HEADER extern bool OCCMeshFace (const OCCGeometry & geom, Mesh & mesh, FlatArray<int, PointIndex> glob2loc,
const MeshingParameters & mparam, int nr, int projecttype, bool delete_on_failure);
namespace step_utils

6
libsrc/occ/occmeshsurf.cpp
View File

@ -2,16 +2,16 @@
#include <mystdlib.h>
#include <occgeom.hpp>
#include <meshing.hpp>
#include "occgeom.hpp"
#include <GeomLProp_SLProps.hxx>
#include <ShapeAnalysis_Surface.hxx>
#include "occmeshsurf.hpp"
namespace netgen
{
#include "occmeshsurf.hpp"
bool glob_testout(false);

12
libsrc/occ/occmeshsurf.hpp
View File

@ -6,9 +6,15 @@
#include "occgeom.hpp"
#include "mydefs.hpp"
#include <TopoDS_Face.hxx>
#include <Geom_Surface.hxx>
#include <ShapeAnalysis.hxx>
#define PARAMETERSPACE -1
#define PLANESPACE 1
namespace netgen
{
class OCCGeometry;
class SingularMatrixException
@ -144,8 +150,8 @@ protected:
class OCCGeometry;
#endif
} // namespace netgen
#endif
#endif

3
libsrc/occ/occpkg.cpp
View File

@ -15,6 +15,9 @@
#include "vsocc.hpp"
#include <TopoDS_Edge.hxx>
#include <IGESControl_Writer.hxx>
// __declspec(dllimport) void AutoColourBcProps(Mesh & mesh, const char *bccolourfile);
// __declspec(dllimport) void GetFaceColours(Mesh & mesh, NgArray<Vec3d> & face_colours);
// __declspec(dllimport) bool ColourMatch(Vec3d col1, Vec3d col2, double eps = 2.5e-05);

61
libsrc/occ/python_occ.cpp
View File

@ -1,58 +1,25 @@
#ifdef NG_PYTHON
#ifdef OCCGEOMETRY
#include <../general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include "../meshing/python_mesh.hpp"
#include <memory>
#include <general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include <meshing/python_mesh.hpp>
#include <meshing.hpp>
#include <occgeom.hpp>
#include <gp_Ax1.hxx>
#include <gp_Ax2.hxx>
#include <gp_Ax2d.hxx>
#include <gp_Trsf.hxx>
#include <BRepPrimAPI_MakeSphere.hxx>
#include <BRepPrimAPI_MakeCylinder.hxx>
#include <BRepPrimAPI_MakeRevol.hxx>
#include <BRepPrimAPI_MakeBox.hxx>
#include <BRepPrimAPI_MakePrism.hxx>
#include <BRepOffsetAPI_MakePipe.hxx>
#include <BRepAlgoAPI_Cut.hxx>
#include <BRepAlgoAPI_Common.hxx>
#include <BRepAlgoAPI_Fuse.hxx>
// #include <XCAFDoc_VisMaterialTool.hxx>
#include <TDF_Attribute.hxx>
#include <Standard_GUID.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <Geom_Plane.hxx>
#include <GC_MakeSegment.hxx>
#include <GC_MakeCircle.hxx>
#include <GC_MakeArcOfCircle.hxx>
#include <GC_MakePlane.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_Transform.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepFilletAPI_MakeFillet.hxx>
#include <BRepOffsetAPI_ThruSections.hxx>
#include <BRepGProp.hxx>
#include <BRepOffsetAPI_MakeThickSolid.hxx>
#include <BRepLib.hxx>
#include <Geom2d_Curve.hxx>
#include <Geom2d_Ellipse.hxx>
#include <Geom2d_TrimmedCurve.hxx>
#include <GCE2d_MakeSegment.hxx>
#include <GCE2d_MakeCircle.hxx>
#include "occgeom.hpp"
#include <BOPAlgo_Builder.hxx>
#include <BRepLProp_SLProps.hxx>
#include <Message.hxx>
#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
#define OCC_HAVE_DUMP_JSON
#endif
#include <Standard_GUID.hxx>
#include <Standard_Version.hxx>
#include <TDF_Attribute.hxx>
#include <XCAFApp_Application.hxx>
#include <XCAFDoc_DocumentTool.hxx>
#include <XCAFDoc_MaterialTool.hxx>
#include <XCAFDoc_ShapeTool.hxx>
using namespace netgen;
@ -297,8 +264,6 @@ DLL_HEADER void ExportNgOCC(py::module &m)
auto result = geo->GenerateMesh(mesh, mp);
if(result != 0)
throw Exception("Meshing failed!");
if (geo->occdim==2)
mesh->SetDimension(2);
SetGlobalMesh(mesh);
ng_geometry = geo;
return mesh;

68
libsrc/occ/python_occ.hpp
View File

@ -1,68 +0,0 @@
class DirectionalInterval
{
public:
gp_Vec dir;
double minval = -1e99;
double maxval = 1e99;
bool openmin = false, openmax = false;
DirectionalInterval (gp_Vec adir) : dir(adir) { ; }
DirectionalInterval (const DirectionalInterval & i2)
: dir(i2.dir), minval(i2.minval), maxval(i2.maxval) { ; }
DirectionalInterval operator< (double val) const
{
DirectionalInterval i2 = *this;
i2.maxval = val;
return i2;
}
DirectionalInterval operator> (double val) const
{
DirectionalInterval i2 = *this;
i2.minval = val;
return i2;
}
DirectionalInterval Intersect (const DirectionalInterval & i2)
{
DirectionalInterval res = *this;
res.minval = max(res.minval, i2.minval);
res.maxval = min(res.maxval, i2.maxval);
return res;
}
bool Contains (gp_Pnt p, double eps = 1e-8)
{
// cout << "Contains point " << p.X() << "," << p.Y() << "," << p.Z() << " ? " << endl;
double val = dir.X()*p.X() + dir.Y()*p.Y() + dir.Z() * p.Z();
// cout << "minval = " << minval << ", val = " << val << " maxval = " << maxval << endl;
if (openmin) {
if (val < minval+eps) return false;
} else {
if (val < minval-eps) return false;
}
if (openmax) {
if (val > maxval-eps) return false;
} else {
if (val > maxval+eps) return false;
}
return true;
}
};
inline gp_Pnt Center (TopoDS_Shape shape)
{
GProp_GProps props;
switch (shape.ShapeType())
{
case TopAbs_FACE:
BRepGProp::SurfaceProperties (shape, props); break;
default:
BRepGProp::LinearProperties(shape, props);
}
return props.CentreOfMass();
}

52
libsrc/occ/python_occ_basic.cpp
View File

@ -1,59 +1,21 @@
#ifdef NG_PYTHON
#ifdef OCCGEOMETRY
#include <../general/ngpython.hpp>
#include <general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include "../meshing/python_mesh.hpp"
#include <meshing/python_mesh.hpp>
#include <meshing.hpp>
#include <occgeom.hpp>
#include "occgeom.hpp"
#include <BRepBuilderAPI_Transform.hxx>
#include <gp_Ax1.hxx>
#include <gp_Ax2.hxx>
#include <gp_Ax2d.hxx>
#include <gp_Ax3.hxx>
#include <gp_Trsf.hxx>
#include <BRepPrimAPI_MakeSphere.hxx>
#include <BRepPrimAPI_MakeCylinder.hxx>
#include <BRepPrimAPI_MakeBox.hxx>
#include <BRepPrimAPI_MakePrism.hxx>
#include <BRepOffsetAPI_MakePipe.hxx>
#include <BRepAlgoAPI_Cut.hxx>
#include <BRepAlgoAPI_Common.hxx>
#include <BRepAlgoAPI_Fuse.hxx>
// #include <XCAFDoc_VisMaterialTool.hxx>
#include <TDF_Attribute.hxx>
#include <Standard_GUID.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <GC_MakeSegment.hxx>
#include <GC_MakeCircle.hxx>
#include <GC_MakeArcOfCircle.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_Transform.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepFilletAPI_MakeFillet.hxx>
#include <BRepOffsetAPI_ThruSections.hxx>
#include <BRepGProp.hxx>
#include <BRepOffsetAPI_MakeThickSolid.hxx>
#include <BRepLib.hxx>
#include <Geom2d_Curve.hxx>
#include <Geom2d_Ellipse.hxx>
#include <Geom2d_TrimmedCurve.hxx>
#include <GCE2d_MakeSegment.hxx>
#include <GCE2d_MakeCircle.hxx>
#include <python_occ.hpp>
#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
#define OCC_HAVE_DUMP_JSON
#endif
using namespace netgen;
DLL_HEADER void ExportNgOCCBasic(py::module &m)
{

284
libsrc/occ/python_occ_shapes.cpp
View File

@ -3,179 +3,68 @@
#include <regex>
#include <../general/ngpython.hpp>
#include <general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include "../meshing/python_mesh.hpp"
#include <meshing/python_mesh.hpp>
#include <meshing.hpp>
#include <occgeom.hpp>
#include <gp_Ax1.hxx>
#include <gp_Ax2.hxx>
#include <gp_Ax2d.hxx>
#include <gp_Trsf.hxx>
#include <BRepPrimAPI_MakeSphere.hxx>
#include <BRepPrimAPI_MakeCylinder.hxx>
#include <BRepPrimAPI_MakeRevol.hxx>
#include <BRepPrimAPI_MakeBox.hxx>
#include <BRepPrimAPI_MakePrism.hxx>
#include <BRepPrimAPI_MakeHalfSpace.hxx>
#include <BRepOffsetAPI_MakePipe.hxx>
#include <BRepOffsetAPI_MakePipeShell.hxx>
#include <BRepAlgoAPI_Cut.hxx>
#include "occgeom.hpp"
#include <BOPAlgo_Builder.hxx>
#include <BOPTools_AlgoTools.hxx>
#include <BRepAlgoAPI_Common.hxx>
#include <BRepAlgoAPI_Cut.hxx>
#include <BRepAlgoAPI_Fuse.hxx>
// #include <XCAFDoc_VisMaterialTool.hxx>
#include <TDF_Attribute.hxx>
#include <TDataStd_Real.hxx>
#include <TDataStd_Name.hxx>
#include <Standard_GUID.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <Geom_Plane.hxx>
#include <Geom_BSplineCurve.hxx>
#include <Geom_BezierCurve.hxx>
#include <GeomAPI_PointsToBSpline.hxx>
#include <GC_MakeSegment.hxx>
#include <GC_MakeCircle.hxx>
#include <GC_MakeArcOfCircle.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepBuilderAPI_MakeVertex.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_Transform.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepFilletAPI_MakeFillet.hxx>
#include <BRepFilletAPI_MakeChamfer.hxx>
#include <BRepOffsetAPI_ThruSections.hxx>
#include <BRepOffsetAPI_MakeOffset.hxx>
#include <BRepExtrema_DistShapeShape.hxx>
#include <BRepFilletAPI_MakeFillet.hxx>
#include <BRepGProp.hxx>
#include <BRepOffsetAPI_MakeThickSolid.hxx>
#include <BRepLProp_SLProps.hxx>
#include <BRepLib.hxx>
#include <BRepMesh_IncrementalMesh.hxx>
#include <BRepOffsetAPI_MakeOffset.hxx>
#include <BRepOffsetAPI_MakePipe.hxx>
#include <BRepOffsetAPI_MakePipeShell.hxx>
#include <BRepOffsetAPI_MakeThickSolid.hxx>
#include <BRepOffsetAPI_ThruSections.hxx>
#include <BRepPrimAPI_MakeBox.hxx>
#include <BRepPrimAPI_MakeCylinder.hxx>
#include <BRepPrimAPI_MakeHalfSpace.hxx>
#include <BRepPrimAPI_MakePrism.hxx>
#include <BRepPrimAPI_MakeRevol.hxx>
#include <BRepPrimAPI_MakeSphere.hxx>
#include <BRepTools.hxx>
#include <GCE2d_MakeArcOfCircle.hxx>
#include <GCE2d_MakeCircle.hxx>
#include <GCE2d_MakeSegment.hxx>
#include <GC_MakeArcOfCircle.hxx>
#include <GC_MakeCircle.hxx>
#include <GC_MakeSegment.hxx>
#include <GProp_GProps.hxx>
#include <Geom2d_Curve.hxx>
#include <Geom2d_Ellipse.hxx>
#include <Geom2d_TrimmedCurve.hxx>
#include <GCE2d_MakeSegment.hxx>
#include <GCE2d_MakeCircle.hxx>
#include <GCE2d_MakeArcOfCircle.hxx>
#include <ShapeUpgrade_UnifySameDomain.hxx>
#include <GeomAPI_PointsToBSpline.hxx>
#include <GeomLProp_SLProps.hxx>
#include <BOPTools_AlgoTools.hxx>
#include <Geom_BSplineCurve.hxx>
#include <Geom_BezierCurve.hxx>
#include <Geom_Plane.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <IntTools_Context.hxx>
#include <STEPControl_Writer.hxx>
#include <ShapeAnalysis_FreeBounds.hxx>
#include <python_occ.hpp>
#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
#define OCC_HAVE_DUMP_JSON
#endif
#include <ShapeUpgrade_UnifySameDomain.hxx>
#include <gp_Ax1.hxx>
#include <gp_Ax2.hxx>
#include <gp_Ax2d.hxx>
#include <gp_Pln.hxx>
#include <gp_Trsf.hxx>
using namespace netgen;
struct ShapeLess
{
bool operator() (const TopoDS_Shape& s1, const TopoDS_Shape& s2) const
{
return s1.TShape() < s2.TShape();
}
};
class ListOfShapes : public std::vector<TopoDS_Shape>
{
public:
TopoDS_Shape Max(gp_Vec dir)
{
double maxval = -1e99;
TopoDS_Shape maxshape;
for (auto shape : *this)
{
GProp_GProps props;
gp_Pnt center;
switch (shape.ShapeType())
{
case TopAbs_VERTEX:
center = BRep_Tool::Pnt (TopoDS::Vertex(shape)); break;
case TopAbs_FACE:
BRepGProp::SurfaceProperties (shape, props);
center = props.CentreOfMass();
break;
default:
BRepGProp::LinearProperties(shape, props);
center = props.CentreOfMass();
}
double val = center.X()*dir.X() + center.Y()*dir.Y() + center.Z() * dir.Z();
if (val > maxval)
{
maxval = val;
maxshape = shape;
}
}
return maxshape;
}
TopoDS_Shape Nearest(gp_Pnt pnt)
{
double mindist = 1e99;
TopoDS_Shape nearestshape;
auto vertex = BRepBuilderAPI_MakeVertex (pnt).Vertex();
for (auto shape : *this)
{
double dist = BRepExtrema_DistShapeShape(shape, vertex).Value();
if (dist < mindist)
{
nearestshape = shape;
mindist = dist;
}
}
return nearestshape;
}
ListOfShapes SubShapes(TopAbs_ShapeEnum type) const
{
std::set<TopoDS_Shape, ShapeLess> unique_shapes;
for(const auto& shape : *this)
for(TopExp_Explorer e(shape, type); e.More(); e.Next())
unique_shapes.insert(e.Current());
ListOfShapes sub;
for(const auto& shape : unique_shapes)
sub.push_back(shape);
return sub;
}
ListOfShapes Solids() const
{
return SubShapes(TopAbs_SOLID);
}
ListOfShapes Faces() const
{
return SubShapes(TopAbs_FACE);
}
ListOfShapes Edges() const
{
return SubShapes(TopAbs_EDGE);
}
ListOfShapes Vertices() const
{
return SubShapes(TopAbs_VERTEX);
}
ListOfShapes operator*(const ListOfShapes& other) const
{
ListOfShapes common;
for(const auto& shape : (*this))
for(const auto& shape_o : other)
if(shape.IsSame(shape_o))
common.push_back(shape);
return common;
}
};
void ExtractEdgeData( const TopoDS_Edge & edge, int index, std::vector<double> * p, Box<3> & box )
{
if (BRep_Tool::Degenerated(edge)) return;
@ -717,36 +606,14 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
return sub;
}, py::arg("type"), "returns list of sub-shapes of type 'type'")
.def_property_readonly("solids", [] (const TopoDS_Shape & shape)
{
ListOfShapes solids;
for(TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
solids.push_back(e.Current());
return solids;
}, "returns all sub-shapes of type 'SOLID'")
.def_property_readonly("faces", [] (const TopoDS_Shape & shape)
{
ListOfShapes sub;
for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
sub.push_back(e.Current());
return sub;
}, "returns all sub-shapes of type 'FACE'")
.def_property_readonly("edges", [] (const TopoDS_Shape & shape)
{
ListOfShapes sub;
for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next())
sub.push_back(e.Current());
return sub;
}, "returns all sub-shapes of type 'EDGE'")
.def_property_readonly("vertices", [] (const TopoDS_Shape & shape)
{
ListOfShapes sub;
for (TopExp_Explorer e(shape, TopAbs_VERTEX); e.More(); e.Next())
sub.push_back(e.Current());
return sub;
}, "returns all sub-shapes of type 'VERTEX'")
.def_property_readonly("solids", GetSolids,
"returns all sub-shapes of type 'SOLID'")
.def_property_readonly("faces", GetFaces,
"returns all sub-shapes of type 'FACE'")
.def_property_readonly("edges", GetEdges,
"returns all sub-shapes of type 'EDGE'")
.def_property_readonly("vertices", GetVertices,
"returns all sub-shapes of type 'VERTEX'")
.def("Properties", [] (const TopoDS_Shape & shape)
{
@ -1127,30 +994,27 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
BRepMesh_IncrementalMesh (shape, deflection, true);
})
.def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name) {
.def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE idtype) {
// only edges supported, by now
auto me_edge = TopoDS::Edge(me);
auto you_edge = TopoDS::Edge(you);
auto type = me.ShapeType();
auto tyou = you.ShapeType();
if(type != tyou)
throw NgException ("Identify: cannot identify different shape types");
GProp_GProps props;
BRepGProp::LinearProperties(me, props);
gp_Pnt cme = props.CentreOfMass();
BRepGProp::LinearProperties(you, props);
gp_Pnt cyou = props.CentreOfMass();
double s0, s1;
auto curve_me = BRep_Tool::Curve(me_edge, s0, s1);
auto vme = occ2ng(curve_me->Value(s1))-occ2ng(curve_me->Value(s0));
auto curve_you = BRep_Tool::Curve(you_edge, s0, s1);
auto vyou = occ2ng(curve_you->Value(s1))-occ2ng(curve_you->Value(s0));
bool inv = vme*vyou < 0;
OCCGeometry::identifications[me.TShape()].push_back
(OCCIdentification { you, Transformation<3>(occ2ng(cyou) - occ2ng(cme)), inv, name });
OCCGeometry::identifications[you.TShape()].push_back
(OCCIdentification { me, Transformation<3>(occ2ng(cme) - occ2ng(cyou)), inv, name });
}, py::arg("other"), py::arg("name"), "Identify shapes for periodic meshing")
switch(type)
{
case TopAbs_VERTEX:
case TopAbs_EDGE:
OCCGeometry::IdentifyEdges(me, you, name, idtype);
break;
default:
throw NgException ("Identify: unsupported shape type");
break;
}
}, py::arg("other"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, "Identify shapes for periodic meshing")
.def("Identify", OCCGeometry::IdentifyFaces, "Identify faces",
py::arg("from"), py::arg("to"), py::arg("name"), py::arg("type")=Identifications::PERIODIC)
.def("Triangulation", [](const TopoDS_Shape & shape)
{
@ -1665,6 +1529,18 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
OCCGeometry::global_shape_properties[shape.TShape()].maxh = maxh;
}
}, "set maxh for all elements of list")
.def_property("hpref", [](ListOfShapes& shapes)
{
throw Exception("Cannot get property of ListOfShapes, get the property from individual shapes!");
},
[](ListOfShapes& shapes, double hpref)
{
for(auto& shape : shapes)
{
auto& val = OCCGeometry::global_shape_properties[shape.TShape()].hpref;
val = max2(hpref, val);
}
}, "set hpref for all elements of list")
;

25
libsrc/occ/vsocc.cpp
View File

@ -8,19 +8,18 @@
#include <occgeom.hpp>
#include "TopoDS_Shape.hxx"
#include "TopoDS_Vertex.hxx"
#include "TopExp_Explorer.hxx"
#include "BRep_Tool.hxx"
#include "TopoDS.hxx"
#include "gp_Pnt.hxx"
#include "Geom_Curve.hxx"
#include "Poly_Triangulation.hxx"
#include "Poly_Array1OfTriangle.hxx"
#include "TColgp_Array1OfPnt2d.hxx"
#include "Poly_Triangle.hxx"
#include "Poly_Polygon3D.hxx"
#include "Poly_PolygonOnTriangulation.hxx"
#include <BRepAdaptor_Surface.hxx>
#include <BRepBndLib.hxx>
#include <BRepLProp_SLProps.hxx>
#include <BRep_Tool.hxx>
#include <Bnd_Box.hxx>
#include <Geom_Curve.hxx>
#include <Poly_PolygonOnTriangulation.hxx>
#include <Poly_Triangle.hxx>
#include <Poly_Triangulation.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Edge.hxx>
#include <gp_Pnt.hxx>
#include <visual.hpp>

6
nglib/nglib.cpp
View File

@ -880,7 +880,7 @@ namespace nglib
mp->Transfer_Parameters();
OCCFindEdges(*occgeom, *me, mparam);
occgeom->FindEdges(*me, mparam);
if((me->GetNP()) && (me->GetNFD()))
{
@ -926,8 +926,8 @@ namespace nglib
perfstepsend = MESHCONST_OPTSURFACE;
}
OCCMeshSurface(*occgeom, *me, mparam);
OCCOptimizeSurface(*occgeom, *me, mparam);
occgeom->MeshSurface(*me, mparam);
occgeom->OptimizeSurface(*me, mparam);
me->CalcSurfacesOfNode();

4
python/occ.py
View File

@ -9,6 +9,10 @@ For more detailed documentation consult the OCCT docs, a good starting point is
https://dev.opencascade.org/doc/refman/html/class_b_rep_builder_a_p_i___make_shape.html
"""
from .config import USE_OCC
if not USE_OCC:
raise ImportError("Netgen was not built with Opencascade support")
from .libngpy._NgOCC import *
from .meshing import meshsize

3
tests/dockerfile
View File

@ -11,12 +11,15 @@ RUN apt-get update && apt-get -y install \
libcgns-dev \
libglu1-mesa-dev \
libhdf5-dev \
libocct-ocaf-dev \
libocct-visualization-dev \
libocct-data-exchange-dev \
libocct-draw-dev \
libpython3-dev \
libtbb-dev \
libxi-dev \
libxmu-dev \
occt-misc \
python3 \
python3-distutils \
python3-numpy \

23
tests/dockerfile_mpi
View File

@ -1,6 +1,27 @@
FROM ubuntu:20.04
ENV DEBIAN_FRONTEND=noninteractive
MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
RUN apt-get update && apt-get -y install python3 libpython3-dev python3-pip libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-numpy python3-tk python3-mpi4py clang-tidy python3-distutils clang libopenmpi-dev openmpi-bin gfortran
RUN apt-get update && apt-get -y install \
ccache \
clang \
clang-tidy \
cmake \
g++ \
gfortran \
git \
libglu1-mesa-dev \
libopenmpi-dev \
libpython3-dev \
libxmu-dev \
openmpi-bin \
python3 \
python3-distutils \
python3-mpi4py \
python3-numpy \
python3-pip \
python3-tk \
tcl-dev \
tk-dev
RUN python3 -m pip install pytest-mpi pytest-check pytest
ADD . /root/src/netgen