mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-25 21:40:33 +05:00
87e74276f3
* Removes dependency on OCC for mesh visualization even for meshes created from OCC Geometry * Code cleanup to remove old commented out functions
1578 lines
42 KiB
C++
1578 lines
42 KiB
C++
#ifdef OCCGEOMETRY
|
|
|
|
#include <mystdlib.h>
|
|
#include <occgeom.hpp>
|
|
#include <meshing.hpp>
|
|
|
|
#include <stlgeom.hpp>
|
|
|
|
|
|
namespace netgen
|
|
{
|
|
|
|
#include "occmeshsurf.hpp"
|
|
|
|
#define TCL_OK 0
|
|
#define TCL_ERROR 1
|
|
|
|
extern STLParameters stlparam;
|
|
|
|
#define DIVIDEEDGESECTIONS 1000
|
|
#define IGNORECURVELENGTH 1e-4
|
|
|
|
|
|
bool merge_solids = 1;
|
|
|
|
inline Point<3> occ2ng (const gp_Pnt & p)
|
|
{
|
|
return Point<3> (p.X(), p.Y(), p.Z());
|
|
}
|
|
|
|
void DivideEdge (TopoDS_Edge & edge,
|
|
Array<MeshPoint> & ps,
|
|
Array<double> & params,
|
|
Mesh & mesh)
|
|
{
|
|
double s0, s1;
|
|
double maxh = mparam.maxh;
|
|
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);
|
|
|
|
double olddist = 0;
|
|
double dist = 0;
|
|
|
|
for (int i = 1; i <= DIVIDEEDGESECTIONS; 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;
|
|
|
|
|
|
olddist = dist;
|
|
dist = pnt.Distance(oldpnt);
|
|
}
|
|
|
|
// 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;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
static void FindEdges (OCCGeometry & geom, Mesh & mesh)
|
|
{
|
|
const char * savetask = multithread.task;
|
|
multithread.task = "Edge meshing";
|
|
|
|
(*testout) << "edge meshing" << endl;
|
|
|
|
int nvertices = geom.vmap.Extent();
|
|
int nedges = geom.emap.Extent();
|
|
|
|
(*testout) << "nvertices = " << nvertices << endl;
|
|
(*testout) << "nedges = " << nedges << endl;
|
|
|
|
double eps = 1e-6 * geom.GetBoundingBox().Diam();
|
|
|
|
for (int i = 1; i <= nvertices; i++)
|
|
{
|
|
gp_Pnt pnt = BRep_Tool::Pnt (TopoDS::Vertex(geom.vmap(i)));
|
|
MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) );
|
|
|
|
bool exists = 0;
|
|
if (merge_solids)
|
|
for (PointIndex pi = 1; pi <= mesh.GetNP(); pi++)
|
|
if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps)
|
|
{
|
|
exists = 1;
|
|
break;
|
|
}
|
|
|
|
if (!exists)
|
|
mesh.AddPoint (mp);
|
|
}
|
|
|
|
(*testout) << "different vertices = " << mesh.GetNP() << endl;
|
|
|
|
|
|
int first_ep = mesh.GetNP()+1;
|
|
|
|
Array<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 (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 facenr = 0;
|
|
int edgenr = 0;
|
|
|
|
|
|
(*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));
|
|
facenr = geom.fmap.FindIndex (face); // sollte doch immer == i3 sein ??? JS
|
|
|
|
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));
|
|
|
|
// Philippose - 06/07/2009
|
|
// Add the face colour to the mesh data
|
|
Quantity_Color face_colour;
|
|
|
|
if(!(geom.face_colours.IsNull())
|
|
&& (geom.face_colours->GetColor(face,XCAFDoc_ColorSurf,face_colour)))
|
|
{
|
|
mesh.GetFaceDescriptor(facenr).SetSurfColour(Vec3d(face_colour.Red(),face_colour.Green(),face_colour.Blue()));
|
|
}
|
|
else
|
|
{
|
|
mesh.GetFaceDescriptor(facenr).SetSurfColour(Vec3d(0.0,1.0,0.0));
|
|
}
|
|
// ACHTUNG! STIMMT NICHT ALLGEMEIN (RG)
|
|
|
|
|
|
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 (TopExp_Explorer exp3 (wire, TopAbs_EDGE); exp3.More(); exp3.Next())
|
|
{
|
|
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.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;
|
|
}
|
|
}
|
|
|
|
|
|
Handle(Geom2d_Curve) cof;
|
|
double s0, s1;
|
|
cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
|
|
|
|
int geomedgenr = geom.emap.FindIndex(edge);
|
|
|
|
Array <MeshPoint> mp;
|
|
Array <double> params;
|
|
|
|
DivideEdge (edge, mp, params, mesh);
|
|
|
|
|
|
Array <int> pnums;
|
|
pnums.SetSize (mp.Size()+2);
|
|
|
|
if (!merge_solids)
|
|
{
|
|
pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge));
|
|
pnums[pnums.Size()-1] = geom.vmap.FindIndex (TopExp::LastVertex (edge));
|
|
}
|
|
else
|
|
{
|
|
Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));
|
|
Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));
|
|
|
|
pnums[0] = -1;
|
|
pnums.Last() = -1;
|
|
for (PointIndex pi = 1; pi < first_ep; pi++)
|
|
{
|
|
if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;
|
|
if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;
|
|
}
|
|
}
|
|
|
|
|
|
for (int i = 1; i <= mp.Size(); i++)
|
|
{
|
|
bool exists = 0;
|
|
int j;
|
|
for (j = first_ep; j <= mesh.GetNP(); j++)
|
|
if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
|
|
{
|
|
exists = 1;
|
|
break;
|
|
}
|
|
|
|
if (exists)
|
|
pnums[i] = j;
|
|
else
|
|
{
|
|
mesh.AddPoint (mp[i-1]);
|
|
(*testout) << "add meshpoint " << mp[i-1] << endl;
|
|
pnums[i] = mesh.GetNP();
|
|
}
|
|
}
|
|
(*testout) << "NP = " << mesh.GetNP() << endl;
|
|
|
|
//(*testout) << pnums[pnums.Size()-1] << endl;
|
|
|
|
for (int i = 1; i <= mp.Size()+1; i++)
|
|
{
|
|
edgenr++;
|
|
Segment seg;
|
|
|
|
seg[0] = pnums[i-1];
|
|
seg[1] = pnums[i];
|
|
seg.edgenr = edgenr;
|
|
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;
|
|
|
|
gp_Pnt2d p2d;
|
|
p2d = cof->Value(params[i-1]);
|
|
// if (i == 1) p2d = cof->Value(s0);
|
|
seg.epgeominfo[0].u = p2d.X();
|
|
seg.epgeominfo[0].v = p2d.Y();
|
|
p2d = cof->Value(params[i]);
|
|
// if (i == mp.Size()+1) p2d = cof -> Value(s1);
|
|
seg.epgeominfo[1].u = p2d.X();
|
|
seg.epgeominfo[1].v = p2d.Y();
|
|
|
|
/*
|
|
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;
|
|
}
|
|
|
|
|
|
|
|
|
|
static void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend)
|
|
{
|
|
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();
|
|
|
|
Array<int> glob2loc(noldp);
|
|
|
|
//int projecttype = PARAMETERSPACE;
|
|
|
|
int projecttype = PARAMETERSPACE;
|
|
|
|
int notrys = 1;
|
|
|
|
int surfmesherror = 0;
|
|
|
|
for (k = 1; k <= mesh.GetNFD(); k++)
|
|
{
|
|
if(1==0 && !geom.fvispar[k-1].IsDrawable())
|
|
{
|
|
(*testout) << "ignoring face " << k << endl;
|
|
cout << "ignoring face " << k << endl;
|
|
continue;
|
|
}
|
|
|
|
(*testout) << "mesh face " << k << endl;
|
|
multithread.percent = 100 * k / (mesh.GetNFD()+1e-10);
|
|
geom.facemeshstatus[k-1] = -1;
|
|
|
|
|
|
/*
|
|
if (k != 42)
|
|
{
|
|
cout << "skipped" << endl;
|
|
continue;
|
|
}
|
|
*/
|
|
|
|
|
|
FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
|
|
|
|
int oldnf = mesh.GetNSE();
|
|
|
|
Box<3> bb = geom.GetBoundingBox();
|
|
|
|
// int projecttype = PLANESPACE;
|
|
|
|
Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype);
|
|
|
|
if (meshing.GetProjectionType() == PLANESPACE)
|
|
PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)");
|
|
else
|
|
PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)");
|
|
|
|
if (surfmesherror)
|
|
cout << "Surface meshing error occured before (in " << surfmesherror << " faces)" << endl;
|
|
|
|
// Meshing2OCCSurfaces meshing(f2, bb);
|
|
meshing.SetStartTime (starttime);
|
|
|
|
//(*testout) << "Face " << k << endl << endl;
|
|
|
|
|
|
if (meshing.GetProjectionType() == PLANESPACE)
|
|
{
|
|
int cntp = 0;
|
|
glob2loc = 0;
|
|
for (i = 1; i <= mesh.GetNSeg(); i++)
|
|
{
|
|
Segment & seg = mesh.LineSegment(i);
|
|
if (seg.si == k)
|
|
{
|
|
for (j = 1; j <= 2; j++)
|
|
{
|
|
int pi = (j == 1) ? seg[0] : seg[1];
|
|
if (!glob2loc.Get(pi))
|
|
{
|
|
meshing.AddPoint (mesh.Point(pi), pi);
|
|
cntp++;
|
|
glob2loc.Elem(pi) = cntp;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
for (i = 1; i <= mesh.GetNSeg(); i++)
|
|
{
|
|
Segment & seg = mesh.LineSegment(i);
|
|
if (seg.si == k)
|
|
{
|
|
PointGeomInfo gi0, gi1;
|
|
gi0.trignum = gi1.trignum = k;
|
|
gi0.u = seg.epgeominfo[0].u;
|
|
gi0.v = seg.epgeominfo[0].v;
|
|
gi1.u = seg.epgeominfo[1].u;
|
|
gi1.v = seg.epgeominfo[1].v;
|
|
|
|
meshing.AddBoundaryElement (glob2loc.Get(seg[0]), glob2loc.Get(seg[1]), gi0, gi1);
|
|
//(*testout) << gi0.u << " " << gi0.v << endl;
|
|
//(*testout) << gi1.u << " " << gi1.v << endl;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
int cntp = 0;
|
|
|
|
for (i = 1; i <= mesh.GetNSeg(); i++)
|
|
if (mesh.LineSegment(i).si == k)
|
|
cntp+=2;
|
|
|
|
|
|
Array< PointGeomInfo > gis;
|
|
|
|
gis.SetAllocSize (cntp);
|
|
gis.SetSize (0);
|
|
|
|
for (i = 1; i <= mesh.GetNSeg(); i++)
|
|
{
|
|
Segment & seg = mesh.LineSegment(i);
|
|
if (seg.si == k)
|
|
{
|
|
PointGeomInfo gi0, gi1;
|
|
gi0.trignum = gi1.trignum = k;
|
|
gi0.u = seg.epgeominfo[0].u;
|
|
gi0.v = seg.epgeominfo[0].v;
|
|
gi1.u = seg.epgeominfo[1].u;
|
|
gi1.v = seg.epgeominfo[1].v;
|
|
|
|
int locpnum[2] = {0, 0};
|
|
|
|
for (j = 0; j < 2; j++)
|
|
{
|
|
PointGeomInfo gi = (j == 0) ? gi0 : gi1;
|
|
|
|
int l;
|
|
for (l = 0; l < gis.Size() && locpnum[j] == 0; l++)
|
|
{
|
|
double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v);
|
|
|
|
if (dist < 1e-10)
|
|
locpnum[j] = l+1;
|
|
}
|
|
|
|
if (locpnum[j] == 0)
|
|
{
|
|
int pi = (j == 0) ? seg[0] : seg[1];
|
|
meshing.AddPoint (mesh.Point(pi), pi);
|
|
|
|
gis.SetSize (gis.Size()+1);
|
|
gis[l] = gi;
|
|
locpnum[j] = l+1;
|
|
}
|
|
}
|
|
|
|
meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1);
|
|
//(*testout) << gi0.u << " " << gi0.v << endl;
|
|
//(*testout) << gi1.u << " " << gi1.v << endl;
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
// Philippose - 15/01/2009
|
|
double maxh = geom.face_maxh[k-1];
|
|
//double maxh = mparam.maxh;
|
|
mparam.checkoverlap = 0;
|
|
// int noldpoints = mesh->GetNP();
|
|
int noldsurfel = mesh.GetNSE();
|
|
|
|
GProp_GProps sprops;
|
|
BRepGProp::SurfaceProperties(TopoDS::Face(geom.fmap(k)),sprops);
|
|
meshing.SetMaxArea(2.*sprops.Mass());
|
|
|
|
MESHING2_RESULT res;
|
|
|
|
try {
|
|
res = meshing.GenerateMesh (mesh, maxh, k);
|
|
}
|
|
|
|
catch (SingularMatrixException)
|
|
{
|
|
(*myerr) << "Singular Matrix" << endl;
|
|
res = MESHING2_GIVEUP;
|
|
}
|
|
|
|
catch (UVBoundsException)
|
|
{
|
|
(*myerr) << "UV bounds exceeded" << endl;
|
|
res = MESHING2_GIVEUP;
|
|
}
|
|
|
|
projecttype = PARAMETERSPACE;
|
|
|
|
if (res != MESHING2_OK)
|
|
{
|
|
if (notrys == 1)
|
|
{
|
|
for (int i = noldsurfel+1; i <= mesh.GetNSE(); i++)
|
|
mesh.DeleteSurfaceElement (i);
|
|
|
|
mesh.Compress();
|
|
|
|
cout << "retry Surface " << k << endl;
|
|
|
|
k--;
|
|
projecttype*=-1;
|
|
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 (i = oldnf+1; i <= mesh.GetNSE(); i++)
|
|
mesh.SurfaceElement(i).SetIndex (k);
|
|
|
|
}
|
|
|
|
ofstream problemfile("occmesh.rep");
|
|
|
|
problemfile << "SURFACEMESHING" << endl << endl;
|
|
|
|
if (surfmesherror)
|
|
{
|
|
cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl;
|
|
cout << "SURFACE MESHING ERROR OCCURED 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();
|
|
}
|
|
|
|
|
|
|
|
|
|
if (multithread.terminate || perfstepsend < MESHCONST_OPTSURFACE)
|
|
return;
|
|
|
|
multithread.task = "Optimizing surface";
|
|
|
|
static int timer_opt2d = NgProfiler::CreateTimer ("Optimization 2D");
|
|
NgProfiler::StartTimer (timer_opt2d);
|
|
|
|
for (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()+1e-10);
|
|
|
|
FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
|
|
|
|
PrintMessage (1, "Optimize Surface ", k);
|
|
for (i = 1; i <= mparam.optsteps2d; i++)
|
|
{
|
|
// (*testout) << "optstep " << i << endl;
|
|
if (multithread.terminate) return;
|
|
|
|
{
|
|
MeshOptimize2dOCCSurfaces meshopt(geom);
|
|
meshopt.SetFaceIndex (k);
|
|
meshopt.SetImproveEdges (0);
|
|
meshopt.SetMetricWeight (0.2);
|
|
meshopt.SetWriteStatus (0);
|
|
|
|
// (*testout) << "EdgeSwapping (mesh, (i > mparam.optsteps2d/2))" << endl;
|
|
meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2));
|
|
}
|
|
|
|
if (multithread.terminate) return;
|
|
{
|
|
MeshOptimize2dOCCSurfaces meshopt(geom);
|
|
meshopt.SetFaceIndex (k);
|
|
meshopt.SetImproveEdges (0);
|
|
meshopt.SetMetricWeight (0.2);
|
|
meshopt.SetWriteStatus (0);
|
|
|
|
// (*testout) << "ImproveMesh (mesh)" << endl;
|
|
meshopt.ImproveMesh (mesh);
|
|
}
|
|
|
|
{
|
|
MeshOptimize2dOCCSurfaces meshopt(geom);
|
|
meshopt.SetFaceIndex (k);
|
|
meshopt.SetImproveEdges (0);
|
|
meshopt.SetMetricWeight (0.2);
|
|
meshopt.SetWriteStatus (0);
|
|
|
|
// (*testout) << "CombineImprove (mesh)" << endl;
|
|
meshopt.CombineImprove (mesh);
|
|
}
|
|
|
|
if (multithread.terminate) return;
|
|
{
|
|
MeshOptimize2dOCCSurfaces meshopt(geom);
|
|
meshopt.SetFaceIndex (k);
|
|
meshopt.SetImproveEdges (0);
|
|
meshopt.SetMetricWeight (0.2);
|
|
meshopt.SetWriteStatus (0);
|
|
|
|
// (*testout) << "ImproveMesh (mesh)" << endl;
|
|
meshopt.ImproveMesh (mesh);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
|
|
mesh.CalcSurfacesOfNode();
|
|
mesh.Compress();
|
|
|
|
NgProfiler::StopTimer (timer_opt2d);
|
|
|
|
multithread.task = savetask;
|
|
}
|
|
|
|
double ComputeH (double kappa)
|
|
{
|
|
double hret;
|
|
kappa *= mparam.curvaturesafety;
|
|
|
|
if (mparam.maxh * kappa < 1)
|
|
hret = mparam.maxh;
|
|
else
|
|
hret = 1 / kappa;
|
|
|
|
if (mparam.maxh < hret)
|
|
hret = mparam.maxh;
|
|
|
|
return (hret);
|
|
}
|
|
|
|
|
|
class Line
|
|
{
|
|
public:
|
|
Point<3> p0, p1;
|
|
double Dist (Line l)
|
|
{
|
|
Vec<3> n = p1-p0;
|
|
Vec<3> q = l.p1-l.p0;
|
|
double nq = n*q;
|
|
|
|
Point<3> p = p0 + 0.5*n;
|
|
double lambda = (p-l.p0)*n / nq;
|
|
|
|
if (lambda >= 0 && lambda <= 1)
|
|
{
|
|
double d = (p-l.p0-lambda*q).Length();
|
|
// if (d < 1e-3) d = 1e99;
|
|
return d;
|
|
}
|
|
else
|
|
return 1e99;
|
|
}
|
|
|
|
double Length ()
|
|
{
|
|
return (p1-p0).Length();
|
|
};
|
|
};
|
|
|
|
/*
|
|
void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
|
|
BRepLProp_SLProps * prop, Mesh & mesh, const double maxside, int depth, double h = 0)
|
|
{
|
|
|
|
|
|
gp_Pnt2d parmid;
|
|
|
|
parmid.SetX(0.3*(par0.X()+par1.X()+par2.X()));
|
|
parmid.SetY(0.3*(par0.Y()+par1.Y()+par2.Y()));
|
|
|
|
if (depth == 0)
|
|
{
|
|
double curvature = 0;
|
|
|
|
prop->SetParameters (parmid.X(), parmid.Y());
|
|
if (!prop->IsCurvatureDefined())
|
|
{
|
|
(*testout) << "curvature not defined!" << endl;
|
|
return;
|
|
}
|
|
curvature = max(fabs(prop->MinCurvature()),
|
|
fabs(prop->MaxCurvature()));
|
|
|
|
prop->SetParameters (par0.X(), par0.Y());
|
|
if (!prop->IsCurvatureDefined())
|
|
{
|
|
(*testout) << "curvature not defined!" << endl;
|
|
return;
|
|
}
|
|
curvature = max(curvature,max(fabs(prop->MinCurvature()),
|
|
fabs(prop->MaxCurvature())));
|
|
|
|
prop->SetParameters (par1.X(), par1.Y());
|
|
if (!prop->IsCurvatureDefined())
|
|
{
|
|
(*testout) << "curvature not defined!" << endl;
|
|
return;
|
|
}
|
|
curvature = max(curvature,max(fabs(prop->MinCurvature()),
|
|
fabs(prop->MaxCurvature())));
|
|
|
|
prop->SetParameters (par2.X(), par2.Y());
|
|
if (!prop->IsCurvatureDefined())
|
|
{
|
|
(*testout) << "curvature not defined!" << endl;
|
|
return;
|
|
}
|
|
curvature = max(curvature,max(fabs(prop->MinCurvature()),
|
|
fabs(prop->MaxCurvature())));
|
|
|
|
//(*testout) << "curvature " << curvature << endl;
|
|
|
|
if (curvature < 1e-3)
|
|
{
|
|
//(*testout) << "curvature too small (" << curvature << ")!" << endl;
|
|
return;
|
|
// return war bis 10.2.05 auskommentiert
|
|
}
|
|
|
|
|
|
|
|
h = ComputeH (curvature+1e-10);
|
|
|
|
if(h < 1e-4*maxside)
|
|
return;
|
|
|
|
|
|
if (h > 30) return;
|
|
}
|
|
|
|
if (h < maxside) // && depth < 5)
|
|
{
|
|
//cout << "\r h " << h << flush;
|
|
gp_Pnt2d pm0;
|
|
gp_Pnt2d pm1;
|
|
gp_Pnt2d pm2;
|
|
|
|
//cout << "h " << h << " maxside " << maxside << " depth " << depth << endl;
|
|
//cout << "par0 " << par0.X() << " " << par0.Y()
|
|
//<< " par1 " << par1.X() << " " << par1.Y()
|
|
// << " par2 " << par2.X() << " " << par2.Y()<< endl;
|
|
|
|
|
|
pm0.SetX(0.5*(par1.X()+par2.X())); pm0.SetY(0.5*(par1.Y()+par2.Y()));
|
|
pm1.SetX(0.5*(par0.X()+par2.X())); pm1.SetY(0.5*(par0.Y()+par2.Y()));
|
|
pm2.SetX(0.5*(par1.X()+par0.X())); pm2.SetY(0.5*(par1.Y()+par0.Y()));
|
|
|
|
RestrictHTriangle (pm0, pm1, pm2, prop, mesh, 0.5*maxside, depth+1, h);
|
|
RestrictHTriangle (par0, pm1, pm2, prop, mesh, 0.5*maxside, depth+1, h);
|
|
RestrictHTriangle (par1, pm0, pm2, prop, mesh, 0.5*maxside, depth+1, h);
|
|
RestrictHTriangle (par2, pm1, pm0, prop, mesh, 0.5*maxside, depth+1, h);
|
|
}
|
|
else
|
|
{
|
|
gp_Pnt pnt;
|
|
Point3d p3d;
|
|
|
|
prop->SetParameters (parmid.X(), parmid.Y());
|
|
pnt = prop->Value();
|
|
p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
|
|
mesh.RestrictLocalH (p3d, h);
|
|
|
|
|
|
prop->SetParameters (par0.X(), par0.Y());
|
|
pnt = prop->Value();
|
|
p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
|
|
mesh.RestrictLocalH (p3d, h);
|
|
|
|
prop->SetParameters (par1.X(), par1.Y());
|
|
pnt = prop->Value();
|
|
p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
|
|
mesh.RestrictLocalH (p3d, h);
|
|
|
|
prop->SetParameters (par2.X(), par2.Y());
|
|
pnt = prop->Value();
|
|
p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
|
|
mesh.RestrictLocalH (p3d, h);
|
|
|
|
}
|
|
}
|
|
*/
|
|
|
|
|
|
void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
|
|
BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0)
|
|
{
|
|
int ls = -1;
|
|
|
|
gp_Pnt pnt0,pnt1,pnt2;
|
|
|
|
prop->SetParameters (par0.X(), par0.Y());
|
|
pnt0 = prop->Value();
|
|
|
|
prop->SetParameters (par1.X(), par1.Y());
|
|
pnt1 = prop->Value();
|
|
|
|
prop->SetParameters (par2.X(), par2.Y());
|
|
pnt2 = prop->Value();
|
|
|
|
double aux;
|
|
double maxside = pnt0.Distance(pnt1);
|
|
ls = 2;
|
|
aux = pnt1.Distance(pnt2);
|
|
if(aux > maxside)
|
|
{
|
|
maxside = aux;
|
|
ls = 0;
|
|
}
|
|
aux = pnt2.Distance(pnt0);
|
|
if(aux > maxside)
|
|
{
|
|
maxside = aux;
|
|
ls = 1;
|
|
}
|
|
|
|
|
|
|
|
gp_Pnt2d parmid;
|
|
|
|
parmid.SetX(0.3*(par0.X()+par1.X()+par2.X()));
|
|
parmid.SetY(0.3*(par0.Y()+par1.Y()+par2.Y()));
|
|
|
|
if (depth%3 == 0)
|
|
{
|
|
double curvature = 0;
|
|
|
|
prop->SetParameters (parmid.X(), parmid.Y());
|
|
if (!prop->IsCurvatureDefined())
|
|
{
|
|
(*testout) << "curvature not defined!" << endl;
|
|
return;
|
|
}
|
|
curvature = max(fabs(prop->MinCurvature()),
|
|
fabs(prop->MaxCurvature()));
|
|
|
|
prop->SetParameters (par0.X(), par0.Y());
|
|
if (!prop->IsCurvatureDefined())
|
|
{
|
|
(*testout) << "curvature not defined!" << endl;
|
|
return;
|
|
}
|
|
curvature = max(curvature,max(fabs(prop->MinCurvature()),
|
|
fabs(prop->MaxCurvature())));
|
|
|
|
prop->SetParameters (par1.X(), par1.Y());
|
|
if (!prop->IsCurvatureDefined())
|
|
{
|
|
(*testout) << "curvature not defined!" << endl;
|
|
return;
|
|
}
|
|
curvature = max(curvature,max(fabs(prop->MinCurvature()),
|
|
fabs(prop->MaxCurvature())));
|
|
|
|
prop->SetParameters (par2.X(), par2.Y());
|
|
if (!prop->IsCurvatureDefined())
|
|
{
|
|
(*testout) << "curvature not defined!" << endl;
|
|
return;
|
|
}
|
|
curvature = max(curvature,max(fabs(prop->MinCurvature()),
|
|
fabs(prop->MaxCurvature())));
|
|
|
|
//(*testout) << "curvature " << curvature << endl;
|
|
|
|
if (curvature < 1e-3)
|
|
{
|
|
//(*testout) << "curvature too small (" << curvature << ")!" << endl;
|
|
return;
|
|
// return war bis 10.2.05 auskommentiert
|
|
}
|
|
|
|
|
|
|
|
h = ComputeH (curvature+1e-10);
|
|
|
|
if(h < 1e-4*maxside)
|
|
return;
|
|
|
|
|
|
if (h > 30) return;
|
|
}
|
|
|
|
if (h < maxside && depth < 10)
|
|
{
|
|
//cout << "\r h " << h << flush;
|
|
gp_Pnt2d pm;
|
|
|
|
//cout << "h " << h << " maxside " << maxside << " depth " << depth << endl;
|
|
//cout << "par0 " << par0.X() << " " << par0.Y()
|
|
//<< " par1 " << par1.X() << " " << par1.Y()
|
|
// << " par2 " << par2.X() << " " << par2.Y()<< endl;
|
|
|
|
if(ls == 0)
|
|
{
|
|
pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y()));
|
|
RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h);
|
|
RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h);
|
|
}
|
|
else if(ls == 1)
|
|
{
|
|
pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y()));
|
|
RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h);
|
|
RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h);
|
|
}
|
|
else if(ls == 2)
|
|
{
|
|
pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y()));
|
|
RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h);
|
|
RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h);
|
|
}
|
|
|
|
}
|
|
else
|
|
{
|
|
gp_Pnt pnt;
|
|
Point3d p3d;
|
|
|
|
prop->SetParameters (parmid.X(), parmid.Y());
|
|
pnt = prop->Value();
|
|
p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
|
|
mesh.RestrictLocalH (p3d, h);
|
|
|
|
p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z());
|
|
mesh.RestrictLocalH (p3d, h);
|
|
|
|
p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z());
|
|
mesh.RestrictLocalH (p3d, h);
|
|
|
|
p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z());
|
|
mesh.RestrictLocalH (p3d, h);
|
|
|
|
//(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;
|
|
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
int OCCGenerateMesh (OCCGeometry & geom,
|
|
Mesh *& mesh,
|
|
int perfstepsstart, int perfstepsend,
|
|
char * optstr)
|
|
{
|
|
// int i, j;
|
|
|
|
multithread.percent = 0;
|
|
|
|
if (perfstepsstart <= MESHCONST_ANALYSE)
|
|
{
|
|
delete mesh;
|
|
mesh = new Mesh();
|
|
mesh->geomtype = Mesh::GEOM_OCC;
|
|
mesh->SetGlobalH (mparam.maxh);
|
|
mesh->SetMinimalH (mparam.minh);
|
|
|
|
Array<double> maxhdom;
|
|
maxhdom.SetSize (geom.NrSolids());
|
|
maxhdom = mparam.maxh;
|
|
|
|
mesh->SetMaxHDomain (maxhdom);
|
|
|
|
Box<3> bb = geom.GetBoundingBox();
|
|
bb.Increase (bb.Diam()/10);
|
|
|
|
mesh->SetLocalH (bb.PMin(), bb.PMax(), 0.5);
|
|
|
|
if (mparam.uselocalh)
|
|
{
|
|
|
|
const char * savetask = multithread.task;
|
|
multithread.percent = 0;
|
|
|
|
mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
|
|
|
|
int nedges = geom.emap.Extent();
|
|
|
|
double maxedgelen = 0;
|
|
double minedgelen = 1e99;
|
|
|
|
multithread.task = "Setting local mesh size (elements per edge)";
|
|
|
|
// setting elements per edge
|
|
|
|
for (int i = 1; i <= nedges && !multithread.terminate; i++)
|
|
{
|
|
TopoDS_Edge e = TopoDS::Edge (geom.emap(i));
|
|
multithread.percent = 100 * (i-1)/double(nedges);
|
|
if (BRep_Tool::Degenerated(e)) continue;
|
|
|
|
GProp_GProps system;
|
|
BRepGProp::LinearProperties(e, system);
|
|
double len = system.Mass();
|
|
|
|
if (len < IGNORECURVELENGTH)
|
|
{
|
|
(*testout) << "ignored" << endl;
|
|
continue;
|
|
}
|
|
|
|
double localh = len/mparam.segmentsperedge;
|
|
double s0, s1;
|
|
|
|
// Philippose - 23/01/2009
|
|
// Find all the parent faces of a given edge
|
|
// and limit the mesh size of the edge based on the
|
|
// mesh size limit of the face
|
|
TopTools_IndexedDataMapOfShapeListOfShape edge_face_map;
|
|
edge_face_map.Clear();
|
|
|
|
TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map);
|
|
const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e);
|
|
|
|
TopTools_ListIteratorOfListOfShape parent_face_list;
|
|
|
|
for(parent_face_list.Initialize(parent_faces); parent_face_list.More(); parent_face_list.Next())
|
|
{
|
|
TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value());
|
|
|
|
int face_index = geom.fmap.FindIndex(parent_face);
|
|
|
|
if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]);
|
|
}
|
|
|
|
Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1);
|
|
|
|
maxedgelen = max (maxedgelen, len);
|
|
minedgelen = min (minedgelen, len);
|
|
|
|
// Philippose - 23/01/2009
|
|
// Modified the calculation of maxj, because the
|
|
// method used so far always results in maxj = 2,
|
|
// which causes the localh to be set only at the
|
|
// starting, mid and end of the edge.
|
|
// Old Algorithm:
|
|
// int maxj = 2 * (int) ceil (localh/len);
|
|
int maxj = max((int) ceil(len/localh), 2);
|
|
|
|
for (int j = 0; j <= maxj; j++)
|
|
{
|
|
gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0));
|
|
mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh);
|
|
}
|
|
}
|
|
|
|
multithread.task = "Setting local mesh size (edge curvature)";
|
|
|
|
// setting edge curvature
|
|
|
|
int nsections = 20;
|
|
|
|
for (int i = 1; i <= nedges && !multithread.terminate; i++)
|
|
{
|
|
double maxcur = 0;
|
|
multithread.percent = 100 * (i-1)/double(nedges);
|
|
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
|
|
if (BRep_Tool::Degenerated(edge)) continue;
|
|
double s0, s1;
|
|
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
|
|
BRepAdaptor_Curve brepc(edge);
|
|
BRepLProp_CLProps prop(brepc, 2, 1e-5);
|
|
|
|
for (int j = 1; j <= nsections; j++)
|
|
{
|
|
double s = s0 + j/(double) nsections * (s1-s0);
|
|
prop.SetParameter (s);
|
|
double curvature = prop.Curvature();
|
|
if(curvature> maxcur) maxcur = curvature;
|
|
|
|
if (curvature >= 1e99)
|
|
continue;
|
|
|
|
gp_Pnt pnt = c->Value (s);
|
|
|
|
mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()),
|
|
ComputeH (fabs(curvature)));
|
|
}
|
|
// (*testout) << "edge " << i << " max. curvature: " << maxcur << endl;
|
|
}
|
|
|
|
multithread.task = "Setting local mesh size (face curvature)";
|
|
|
|
// setting face curvature
|
|
|
|
int nfaces = geom.fmap.Extent();
|
|
|
|
for (int i = 1; i <= nfaces && !multithread.terminate; i++)
|
|
{
|
|
multithread.percent = 100 * (i-1)/double(nfaces);
|
|
TopoDS_Face face = TopoDS::Face(geom.fmap(i));
|
|
TopLoc_Location loc;
|
|
Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
|
|
Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
|
|
if (triangulation.IsNull()) continue;
|
|
|
|
BRepAdaptor_Surface sf(face, Standard_True);
|
|
BRepLProp_SLProps prop(sf, 2, 1e-5);
|
|
|
|
int ntriangles = triangulation -> NbTriangles();
|
|
for (int j = 1; j <= ntriangles; j++)
|
|
{
|
|
gp_Pnt p[3];
|
|
gp_Pnt2d par[3];
|
|
|
|
for (int k = 1; k <=3; k++)
|
|
{
|
|
int n = triangulation->Triangles()(j)(k);
|
|
p[k-1] = triangulation->Nodes()(n).Transformed(loc);
|
|
par[k-1] = triangulation->UVNodes()(n);
|
|
}
|
|
|
|
//double maxside = 0;
|
|
//maxside = max (maxside, p[0].Distance(p[1]));
|
|
//maxside = max (maxside, p[0].Distance(p[2]));
|
|
//maxside = max (maxside, p[1].Distance(p[2]));
|
|
//cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;
|
|
|
|
RestrictHTriangle (par[0], par[1], par[2], &prop, *mesh, 0);
|
|
//cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
|
|
}
|
|
}
|
|
|
|
// setting close edges
|
|
|
|
if (stlparam.resthcloseedgeenable)
|
|
{
|
|
multithread.task = "Setting local mesh size (close edges)";
|
|
|
|
int sections = 100;
|
|
|
|
Array<Line> lines(sections*nedges);
|
|
|
|
Box3dTree* searchtree =
|
|
new Box3dTree (bb.PMin(), bb.PMax());
|
|
|
|
int nlines = 0;
|
|
for (int i = 1; i <= nedges && !multithread.terminate; i++)
|
|
{
|
|
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
|
|
if (BRep_Tool::Degenerated(edge)) continue;
|
|
|
|
double s0, s1;
|
|
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
|
|
BRepAdaptor_Curve brepc(edge);
|
|
BRepLProp_CLProps prop(brepc, 1, 1e-5);
|
|
prop.SetParameter (s0);
|
|
|
|
gp_Vec d0 = prop.D1().Normalized();
|
|
double s_start = s0;
|
|
int count = 0;
|
|
for (int j = 1; j <= sections; j++)
|
|
{
|
|
double s = s0 + (s1-s0)*(double)j/(double)sections;
|
|
prop.SetParameter (s);
|
|
gp_Vec d1 = prop.D1().Normalized();
|
|
double cosalpha = fabs(d0*d1);
|
|
if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI)))
|
|
{
|
|
count++;
|
|
gp_Pnt p0 = c->Value (s_start);
|
|
gp_Pnt p1 = c->Value (s);
|
|
lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z());
|
|
lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z());
|
|
|
|
Box3d box;
|
|
box.SetPoint (Point3d(lines[nlines].p0));
|
|
box.AddPoint (Point3d(lines[nlines].p1));
|
|
|
|
searchtree->Insert (box.PMin(), box.PMax(), nlines+1);
|
|
nlines++;
|
|
|
|
s_start = s;
|
|
d0 = d1;
|
|
}
|
|
}
|
|
}
|
|
|
|
Array<int> linenums;
|
|
|
|
for (int i = 0; i < nlines; i++)
|
|
{
|
|
multithread.percent = (100*i)/double(nlines);
|
|
Line & line = lines[i];
|
|
|
|
Box3d box;
|
|
box.SetPoint (Point3d(line.p0));
|
|
box.AddPoint (Point3d(line.p1));
|
|
double maxhline = max (mesh->GetH(box.PMin()),
|
|
mesh->GetH(box.PMax()));
|
|
box.Increase(maxhline);
|
|
|
|
double mindist = 1e99;
|
|
linenums.SetSize(0);
|
|
searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums);
|
|
|
|
for (int j = 0; j < linenums.Size(); j++)
|
|
{
|
|
int num = linenums[j]-1;
|
|
if (i == num) continue;
|
|
if ((line.p0-lines[num].p0).Length2() < 1e-15) continue;
|
|
if ((line.p0-lines[num].p1).Length2() < 1e-15) continue;
|
|
if ((line.p1-lines[num].p0).Length2() < 1e-15) continue;
|
|
if ((line.p1-lines[num].p1).Length2() < 1e-15) continue;
|
|
mindist = min (mindist, line.Dist(lines[num]));
|
|
}
|
|
|
|
mindist *= stlparam.resthcloseedgefac;
|
|
|
|
if (mindist < 1e-3)
|
|
{
|
|
(*testout) << "extremely small local h: " << mindist
|
|
<< " --> setting to 1e-3" << endl;
|
|
(*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl;
|
|
mindist = 1e-3;
|
|
}
|
|
|
|
mesh->RestrictLocalHLine(line.p0, line.p1, mindist);
|
|
}
|
|
}
|
|
|
|
multithread.task = savetask;
|
|
|
|
}
|
|
|
|
// Philippose - 09/03/2009
|
|
// Added the capability to load the mesh size from a
|
|
// file also for OpenCascade Geometry
|
|
// Note:
|
|
// ** If the "uselocalh" option is ticked in
|
|
// the "mesh options...insider" menu, the mesh
|
|
// size will be further modified by the topology
|
|
// analysis routines.
|
|
// ** To use the mesh size file as the sole source
|
|
// for defining the mesh size, uncheck the "uselocalh"
|
|
// option.
|
|
mesh->LoadLocalMeshSize (mparam.meshsizefilename);
|
|
}
|
|
|
|
if (multithread.terminate || perfstepsend <= MESHCONST_ANALYSE)
|
|
return TCL_OK;
|
|
|
|
if (perfstepsstart <= MESHCONST_MESHEDGES)
|
|
{
|
|
FindEdges (geom, *mesh);
|
|
|
|
/*
|
|
cout << "Removing redundant points" << endl;
|
|
|
|
int i, j;
|
|
int np = mesh->GetNP();
|
|
Array<int> equalto;
|
|
|
|
equalto.SetSize (np);
|
|
equalto = 0;
|
|
|
|
for (i = 1; i <= np; i++)
|
|
{
|
|
for (j = i+1; j <= np; j++)
|
|
{
|
|
if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12))
|
|
equalto[j-1] = i;
|
|
}
|
|
}
|
|
|
|
for (i = 1; i <= np; i++)
|
|
if (equalto[i-1])
|
|
{
|
|
cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl;
|
|
for (j = 1; j <= mesh->GetNSeg(); j++)
|
|
{
|
|
Segment & seg = mesh->LineSegment(j);
|
|
if (seg[0] == i) seg[0] = equalto[i-1];
|
|
if (seg[1] == i) seg[1] = equalto[i-1];
|
|
}
|
|
}
|
|
|
|
cout << "Removing degenerated segments" << endl;
|
|
for (j = 1; j <= mesh->GetNSeg(); j++)
|
|
{
|
|
Segment & seg = mesh->LineSegment(j);
|
|
if (seg[0] == seg[1])
|
|
{
|
|
mesh->DeleteSegment(j);
|
|
cout << "Deleting Segment " << j << endl;
|
|
}
|
|
}
|
|
|
|
mesh->Compress();
|
|
*/
|
|
|
|
/*
|
|
for (int i = 1; i <= geom.fmap.Extent(); i++)
|
|
{
|
|
Handle(Geom_Surface) hf1 =
|
|
BRep_Tool::Surface(TopoDS::Face(geom.fmap(i)));
|
|
for (int j = i+1; j <= geom.fmap.Extent(); j++)
|
|
{
|
|
Handle(Geom_Surface) hf2 =
|
|
BRep_Tool::Surface(TopoDS::Face(geom.fmap(j)));
|
|
if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl;
|
|
}
|
|
}
|
|
*/
|
|
|
|
#ifdef LOG_STREAM
|
|
(*logout) << "Edges meshed" << endl
|
|
<< "time = " << GetTime() << " sec" << endl
|
|
<< "points: " << mesh->GetNP() << endl;
|
|
#endif
|
|
}
|
|
|
|
if (multithread.terminate || perfstepsend <= MESHCONST_MESHEDGES)
|
|
return TCL_OK;
|
|
|
|
if (perfstepsstart <= MESHCONST_MESHSURFACE)
|
|
{
|
|
OCCMeshSurface (geom, *mesh, perfstepsend);
|
|
if (multithread.terminate) return TCL_OK;
|
|
|
|
#ifdef LOG_STREAM
|
|
(*logout) << "Surfaces meshed" << endl
|
|
<< "time = " << GetTime() << " sec" << endl
|
|
<< "points: " << mesh->GetNP() << endl;
|
|
#endif
|
|
|
|
#ifdef STAT_STREAM
|
|
(*statout) << mesh->GetNSeg() << " & "
|
|
<< mesh->GetNSE() << " & - &"
|
|
<< GetTime() << " & " << endl;
|
|
#endif
|
|
|
|
// MeshQuality2d (*mesh);
|
|
mesh->CalcSurfacesOfNode();
|
|
}
|
|
|
|
if (multithread.terminate || perfstepsend <= MESHCONST_OPTSURFACE)
|
|
return TCL_OK;
|
|
|
|
if (perfstepsstart <= MESHCONST_MESHVOLUME)
|
|
{
|
|
multithread.task = "Volume meshing";
|
|
|
|
MESHING3_RESULT res =
|
|
MeshVolume (mparam, *mesh);
|
|
|
|
ofstream problemfile("occmesh.rep",ios_base::app);
|
|
|
|
problemfile << "VOLUMEMESHING" << endl << endl;
|
|
if(res != MESHING3_OK)
|
|
problemfile << "ERROR" << endl << endl;
|
|
else
|
|
problemfile << "OK" << endl
|
|
<< mesh->GetNE() << " elements" << endl << endl;
|
|
|
|
problemfile.close();
|
|
|
|
if (res != MESHING3_OK) return TCL_ERROR;
|
|
|
|
if (multithread.terminate) return TCL_OK;
|
|
|
|
RemoveIllegalElements (*mesh);
|
|
if (multithread.terminate) return TCL_OK;
|
|
|
|
MeshQuality3d (*mesh);
|
|
|
|
#ifdef STAT_STREAM
|
|
(*statout) << GetTime() << " & ";
|
|
#endif
|
|
|
|
#ifdef LOG_STREAM
|
|
(*logout) << "Volume meshed" << endl
|
|
<< "time = " << GetTime() << " sec" << endl
|
|
<< "points: " << mesh->GetNP() << endl;
|
|
#endif
|
|
}
|
|
|
|
if (multithread.terminate || perfstepsend <= MESHCONST_MESHVOLUME)
|
|
return TCL_OK;
|
|
|
|
if (perfstepsstart <= MESHCONST_OPTVOLUME)
|
|
{
|
|
multithread.task = "Volume optimization";
|
|
|
|
OptimizeVolume (mparam, *mesh);
|
|
if (multithread.terminate) return TCL_OK;
|
|
|
|
#ifdef STAT_STREAM
|
|
(*statout) << GetTime() << " & "
|
|
<< mesh->GetNE() << " & "
|
|
<< mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl;
|
|
#endif
|
|
|
|
#ifdef LOG_STREAM
|
|
(*logout) << "Volume optimized" << endl
|
|
<< "time = " << GetTime() << " sec" << endl
|
|
<< "points: " << mesh->GetNP() << endl;
|
|
#endif
|
|
|
|
// cout << "Optimization complete" << endl;
|
|
|
|
}
|
|
|
|
(*testout) << "NP: " << mesh->GetNP() << endl;
|
|
for (int i = 1; i <= mesh->GetNP(); i++)
|
|
(*testout) << mesh->Point(i) << endl;
|
|
|
|
(*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl;
|
|
for (int i = 1; i <= mesh->GetNSeg(); i++)
|
|
(*testout) << mesh->LineSegment(i) << endl;
|
|
|
|
return TCL_OK;
|
|
}
|
|
}
|
|
|
|
#endif
|