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

Merge branch 'optimize_occ_createlocalh' into 'master'

Browse Source 
some optimizations for CalcLocalH in occ mesher

See merge request jschoeberl/netgen!267
This commit is contained in:
Joachim Schöberl 2019-10-04 13:05:43 +00:00
commit 48aa301558
1 changed files with 38 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
libsrc/occ/occgenmesh.cpp
View File

@ -69,7 +69,7 @@ namespace netgen
void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h, BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h,
const MeshingParameters & mparam) const MeshingParameters & mparam)
{ {
int ls = -1; int ls = -1;
@ -112,41 +112,41 @@ namespace netgen
{ {
double curvature = 0; double curvature = 0;
prop->SetParameters (parmid.X(), parmid.Y()); prop2->SetParameters (parmid.X(), parmid.Y());
if (!prop->IsCurvatureDefined()) if (!prop2->IsCurvatureDefined())
{ {
(*testout) << "curvature not defined!" << endl; (*testout) << "curvature not defined!" << endl;
return; return;
} }
curvature = max(fabs(prop->MinCurvature()), curvature = max(fabs(prop2->MinCurvature()),
fabs(prop->MaxCurvature())); fabs(prop2->MaxCurvature()));
prop->SetParameters (par0.X(), par0.Y()); prop2->SetParameters (par0.X(), par0.Y());
if (!prop->IsCurvatureDefined()) if (!prop2->IsCurvatureDefined())
{ {
(*testout) << "curvature not defined!" << endl; (*testout) << "curvature not defined!" << endl;
return; return;
} }
curvature = max(curvature,max(fabs(prop->MinCurvature()), curvature = max(curvature,max(fabs(prop2->MinCurvature()),
fabs(prop->MaxCurvature()))); fabs(prop2->MaxCurvature())));
prop->SetParameters (par1.X(), par1.Y()); prop2->SetParameters (par1.X(), par1.Y());
if (!prop->IsCurvatureDefined()) if (!prop2->IsCurvatureDefined())
{ {
(*testout) << "curvature not defined!" << endl; (*testout) << "curvature not defined!" << endl;
return; return;
} }
curvature = max(curvature,max(fabs(prop->MinCurvature()), curvature = max(curvature,max(fabs(prop2->MinCurvature()),
fabs(prop->MaxCurvature()))); fabs(prop2->MaxCurvature())));
prop->SetParameters (par2.X(), par2.Y()); prop2->SetParameters (par2.X(), par2.Y());
if (!prop->IsCurvatureDefined()) if (!prop2->IsCurvatureDefined())
{ {
(*testout) << "curvature not defined!" << endl; (*testout) << "curvature not defined!" << endl;
return; return;
} }
curvature = max(curvature,max(fabs(prop->MinCurvature()), curvature = max(curvature,max(fabs(prop2->MinCurvature()),
fabs(prop->MaxCurvature()))); fabs(prop2->MaxCurvature())));
//(*testout) << "curvature " << curvature << endl; //(*testout) << "curvature " << curvature << endl;
@ -165,7 +165,7 @@ namespace netgen
return; return;
if (h > 30) return; // if (h > 30) return;
} }
if (h < maxside && depth < 10) if (h < maxside && depth < 10)
@ -181,20 +181,20 @@ namespace netgen
if(ls == 0) if(ls == 0)
{ {
pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y())); 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, mparam); RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h, mparam); RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam);
} }
else if(ls == 1) else if(ls == 1)
{ {
pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y())); 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, mparam); RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h, mparam); RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam);
} }
else if(ls == 2) else if(ls == 2)
{ {
pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y())); 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, mparam); RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h, mparam); RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam);
} }
} }
@ -994,6 +994,8 @@ namespace netgen
void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh, void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh,
const MeshingParameters & mparam, const OCCParameters& occparam) const MeshingParameters & mparam, const OCCParameters& occparam)
{ {
static Timer t1("OCCSetLocalMeshSize");
RegionTimer regt(t1);
mesh.SetGlobalH (mparam.maxh); mesh.SetGlobalH (mparam.maxh);
mesh.SetMinimalH (mparam.minh); mesh.SetMinimalH (mparam.minh);
@ -1029,8 +1031,15 @@ namespace netgen
multithread.task = "Setting local mesh size (elements per edge)"; multithread.task = "Setting local mesh size (elements per edge)";
// setting elements per edge // 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);
// setting elements per edge
for (int i = 1; i <= nedges && !multithread.terminate; i++) for (int i = 1; i <= nedges && !multithread.terminate; i++)
{ {
TopoDS_Edge e = TopoDS::Edge (geom.emap(i)); TopoDS_Edge e = TopoDS::Edge (geom.emap(i));
@ -1050,14 +1059,6 @@ namespace netgen
double localh = len/mparam.segmentsperedge; double localh = len/mparam.segmentsperedge;
double s0, s1; 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); const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e);
TopTools_ListIteratorOfListOfShape parent_face_list; TopTools_ListIteratorOfListOfShape parent_face_list;
@ -1140,7 +1141,9 @@ namespace netgen
} }
BRepAdaptor_Surface sf(face, Standard_True); BRepAdaptor_Surface sf(face, Standard_True);
BRepLProp_SLProps prop(sf, 2, 1e-5); // one prop for evaluating and one for derivatives
BRepLProp_SLProps prop(sf, 0, 1e-5);
BRepLProp_SLProps prop2(sf, 2, 1e-5);
int ntriangles = triangulation -> NbTriangles(); int ntriangles = triangulation -> NbTriangles();
for (int j = 1; j <= ntriangles; j++) for (int j = 1; j <= ntriangles; j++)
@ -1161,7 +1164,7 @@ namespace netgen
//maxside = max (maxside, p[1].Distance(p[2])); //maxside = max (maxside, p[1].Distance(p[2]));
//cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush; //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;
RestrictHTriangle (par[0], par[1], par[2], &prop, mesh, 0, 0, mparam); RestrictHTriangle (par[0], par[1], par[2], &prop, &prop2, mesh, 0, 0, mparam);
//cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush; //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
} }
} }