Exception not handled

nicoloc · February 19, 2024, 3:11pm

Hello,
I’m trying to get the mesh of some CAD parts from their step files, using netgen. This works great most of the time, however in some cases the process abruptly ends with not much information: Process finished with exit code 139 (interrupted by signal 11:SIGSEGV).

To understand what’s the problem, I ran more or less the same program in c++. I found out that the problem lies in setting the face curvature in Analyse(*mesh, mparam) → OCCSetLocalMeshSize(*this, mesh, mparam, occparam); In particular, I get an exception Exception("OCC-Triangulation could not be built. Do you have a bounded shape?");

I was wondering if there’s a way of throwing a python exception which could then be handled from python, rather than abruptly ending the program this way (the only way of handling it now in python would be to call GenerateMesh in a separate subprocess and checking for its status once it finishes).

Also, I’ve been able to identify the problematic face. Using the following function in FreeCAD:

face_mesh = MeshPart.meshFromShape(
            Shape=f,
            LinearDeflection=0.01,
            AngularDeflection=0.5,
            Relative=True,
        )

gives indeed an empty mesh for such face. However, after setting Relative=False, I’m able to get a non-empty mesh. I think the reason that the OCC-Triangulation is null in netgen could derive from the fact that BRepMesh_IncrementalMesh (geom.shape, 0.01, true) is being called with isRelative = true (but maybe I’m wrong, I don’t know exactly how these functions work under the hood).

What’s the rationale behind the choice of those parameters in Netgen (i.e., why theLinDeflection = 0.01, and isRelative = true in BRepMesh_IncrementalMesh)? Would it be possible to set those parameters in the MeshingParameters?

christopher · February 19, 2024, 4:49pm

Hi, I added this exception I think last week or 2 ago, because I ran into exactly the same problem, so I assume you were testing in python with an older ngsolve version. The exception should actually be forwarded to python the same way. If you reproduce the behaviour on the newest master in python then there might be a bug, can you test that?

I think most of these parameters were chosen at some early point of development of the occ interface on a range of test cases, but only adapted if needed. Can you provide explanations for them and possible downsides? Preferably would be of course stable parameter but also performance should not suffer.

We can think about an interface on how to make them adaptable from python though.

nicoloc · February 19, 2024, 8:40pm

Thanks for the quick reply!

What version are you referring to? I’m using the python library from a virtual environment where I installed it pip install --upgrade ngsolve. The behavior I reported above was with netgen-mesher 6.2.2307. After updating it to netgen-mesher-6.2.2401, the behavior is still the same (I get the segmentation fault from python). Should I build netgen again from scratch?

Regarding your questions about the parameters in BRepMesh_IncrementalMesh, I think you can find a good explanation on the OpenCascade website. From my understanding, lower values of linear and angular deflection would give you a more accurate tessellation, at the price of increased runtime. I don’t see the parameter isRelative mentioned there, but if you check the documentation it says: if TRUE deflection used for discretization of each edge will be <theLinDeflection> * <size of edge>. Deflection used for the faces will be the maximum deflection of their edges.
I think it would be nice to have control over those parameters (instead of hardcoded ones as it is right now). Probably the most natural place in netgen to set them would be in OCCParameters, and then modify the call to BRepMesh_IncrementalMesh to use them. I’m not sure how those could be managed in python, maybe include them in MeshingParameters. The list is already quite long though and those parameters are only relevant if you’re using OCC.

christopher · February 19, 2024, 10:29pm

this is the relevant commit: throw exception if surface triangulation cannot be built by occ · NGSolve/netgen@ddc50aa · GitHub

It’s not yet included in any release but you can install nightly ngsolve version with

pip install --pre ngsolve

actually I think we only need the tesselation for setting the mesh size. Maybe we can find a method that tries with finer parameters if tesselation building fails. Or occ has something like this somehow built in…

joachim · February 20, 2024, 6:07am

Hi Nicoloc,

thank you for the report and detailed information.

Can we get more information why OCC - Triangulation fails in this case, e.g. from the OCC status code ?

Could you share your example ?

This new exception suggests that the shape is unbounded, what is certainly not the case here (otherwise Relative=false could not generate a triangulation as well).
We should check shape.Infinite and give a definite error message if the shape is infinite.

We could easily introduce one more meshing parameter, but I don’t want to have tons of parameters for which nobody can explain their relevance to Netgen.

Could it be that one of the edges of the face has zero-size, and thus the relative-deflection is set to 0 (according to the documentation) ?

Joachim

nicoloc · February 20, 2024, 9:51am

Just updated it with pip install --pre ngsolve and it works great. Now I get the exception in python:
netgen.libngpy._meshing.NgException: OCC-Triangulation could not be built. Do you have a bounded shape?

nicoloc · February 20, 2024, 11:01am

Hi Joachim,

Unfortunately I’m not allowed to share the step file as it’s under confidential information. It’s a relatively simple CAD part with 6 faces. I tried to gather some information through FreeCAD:

p = Part.read(cad_f)
p.isInfinite()
Out[2]: False
f = p.Faces[3]  # this is the problematic face
f.Area
Out[4]: 16.237798546387225
[e.Length for e in f.Edges]
Out[7]: 
[47.39898578042153,
 3.8007894453146425,
 11.167965277525747,
 8.333213170139574,
 3.0721468452084655,
 7.539822368615504,
 27.28597626138821,
 42.197548223696415,
 16.354753098347285,
 75.73575564851302,
 55.58014254418251,
 3.800860005983231]

I’m also trying to use netgen in c++, but I’m not sure how to get the OCC status code. Can you point me to the right direction? This is roughly the code I’m using:

  auto mesh = std::make_shared<netgen::Mesh>();
  OCCGeometry* geom = LoadOCC_STEP(filename);
  geom->HealGeometry();

  auto sh = geom->GetShape();
  std::cout << "infinite: " << sh.Infinite() << "\n";

  TopoDS_Face face = TopoDS::Face(geom->fmap(4));  // the problematic face
  TopLoc_Location loc;
  Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
  Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);

  if (triangulation.IsNull())
  {
    BRepTools::Clean (geom->shape);
    BRepMesh_IncrementalMesh (geom->shape, 0.01, true);
    triangulation = BRep_Tool::Triangulation (face, loc);
  }

  mesh->SetGeometry(std::shared_ptr<netgen::OCCGeometry>(geom, &NOOP_Deleter));
  auto result = geom->GenerateMesh(mesh, ngParam);

which gives me:

Starting geometry healing procedure (tolerance: 0)
-----------------------------------

- repairing faces
-----------------------------------
Compounds       : 0 (0)
Composite solids: 0 (0)
Solids          : 1 (1)
Shells          : 1 (1)
Wires           : 8 (8)
Faces           : 6 (6)
Edges           : 25 (25)
Vertices        : 19 (19)

Total surface area : 14688.5 (14688.5)

infinite: 0
triangulation entity: 0
terminate called after throwing an instance of 'ngcore::Exception'
  what():  OCC-Triangulation could not be built. Do you have a bounded shape?
[1]    928384 IOT instruction (core dumped)  ./MyMeshingApp

If I change isRelative to false with
BRepMesh_IncrementalMesh (geom->shape, 0.01, false);
I then get:

meshed area = 58.3703
maximal area = 32.4756
GIVING UP
meshed area = 58.3703
maximal area = 32.4756
GIVING UP
 ERROR: Problem in Surface mesh generation
WARNING! NOT ALL FACES HAVE BEEN MESHED
SURFACE MESHING ERROR OCCURRED IN 1 FACES:
 WARNING: Intersecting elements 20 and 1081
el1 = np = 3 3 31 165
el2 = np = 3 19 140 620
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 20 and 1097
el1 = np = 3 3 31 165
el2 = np = 3 19 620 623
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 20 and 1084
el1 = np = 3 3 31 165
el2 = np = 3 144 19 623
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 21 and 1083
el1 = np = 3 31 32 166
el2 = np = 3 143 144 622
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 21 and 1084
el1 = np = 3 31 32 166
el2 = np = 3 144 19 623
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 21 and 1085
el1 = np = 3 31 32 166
el2 = np = 3 144 623 622
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 23 and 1083
el1 = np = 3 32 33 167
el2 = np = 3 143 144 622
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 23 and 1104
el1 = np = 3 32 33 167
el2 = np = 3 142 143 630
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 23 and 1107
el1 = np = 3 32 33 167
el2 = np = 3 630 143 622
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 25 and 157
el1 = np = 3 33 34 168
el2 = np = 3 8 75 224
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 25 and 158
el1 = np = 3 33 34 168
el2 = np = 3 80 8 224
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 25 and 1087
el1 = np = 3 33 34 168
el2 = np = 3 75 8 624
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 25 and 1104
el1 = np = 3 33 34 168
el2 = np = 3 142 143 630
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 25 and 1181
el1 = np = 3 33 34 168
el2 = np = 3 8 142 624
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 25 and 1334
el1 = np = 3 33 34 168
el2 = np = 3 624 142 630
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 27 and 268
el1 = np = 3 34 35 169
el2 = np = 3 79 80 224
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 27 and 158
el1 = np = 3 34 35 169
el2 = np = 3 80 8 224
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 29 and 268
el1 = np = 3 35 36 170
el2 = np = 3 79 80 224
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 29 and 358
el1 = np = 3 35 36 170
el2 = np = 3 79 224 223
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 29 and 163
el1 = np = 3 35 36 170
el2 = np = 3 78 79 227
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 29 and 266
el1 = np = 3 35 36 170
el2 = np = 3 227 79 223
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 31 and 163
el1 = np = 3 36 4 40
el2 = np = 3 78 79 227
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 31 and 211
el1 = np = 3 36 4 40
el2 = np = 3 78 227 251
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 31 and 210
el1 = np = 3 36 4 40
el2 = np = 3 77 78 251
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 157 and 25
el1 = np = 3 8 75 224
el2 = np = 3 33 34 168
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 158 and 27
el1 = np = 3 80 8 224
el2 = np = 3 34 35 169
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 158 and 25
el1 = np = 3 80 8 224
el2 = np = 3 33 34 168
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 161 and 1684
el1 = np = 3 76 77 226
el2 = np = 3 40 4 145
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 161 and 1711
el1 = np = 3 76 77 226
el2 = np = 3 145 9 900
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 162 and 1711
el1 = np = 3 76 226 225
el2 = np = 3 145 9 900
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 163 and 29
el1 = np = 3 78 79 227
el2 = np = 3 35 36 170
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 163 and 31
el1 = np = 3 78 79 227
el2 = np = 3 36 4 40
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 166 and 1715
el1 = np = 3 11 81 229
el2 = np = 3 147 148 903
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 166 and 1713
el1 = np = 3 11 81 229
el2 = np = 3 146 147 902
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 167 and 1713
el1 = np = 3 81 10 230
el2 = np = 3 146 147 902
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 168 and 1717
el1 = np = 3 12 82 99
el2 = np = 3 148 149 904
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 168 and 1719
el1 = np = 3 12 82 99
el2 = np = 3 149 150 905
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 169 and 1717
el1 = np = 3 82 11 231
el2 = np = 3 148 149 904
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 169 and 1715
el1 = np = 3 82 11 231
el2 = np = 3 147 148 903
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 170 and 1717
el1 = np = 3 99 82 231
el2 = np = 3 148 149 904
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 171 and 1715
el1 = np = 3 11 229 231
el2 = np = 3 147 148 903
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 210 and 31
el1 = np = 3 77 78 251
el2 = np = 3 36 4 40
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 210 and 1684
el1 = np = 3 77 78 251
el2 = np = 3 40 4 145
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 211 and 31
el1 = np = 3 78 227 251
el2 = np = 3 36 4 40
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 213 and 1713
el1 = np = 3 229 81 230
el2 = np = 3 146 147 902
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 266 and 29
el1 = np = 3 227 79 223
el2 = np = 3 35 36 170
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 268 and 29
el1 = np = 3 79 80 224
el2 = np = 3 35 36 170
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 268 and 27
el1 = np = 3 79 80 224
el2 = np = 3 34 35 169
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 316 and 1684
el1 = np = 3 77 251 226
el2 = np = 3 40 4 145
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 358 and 29
el1 = np = 3 79 224 223
el2 = np = 3 35 36 170
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1053 and 1719
el1 = np = 3 12 99 17
el2 = np = 3 149 150 905
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1053 and 1721
el1 = np = 3 12 99 17
el2 = np = 3 150 151 906
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1081 and 1676
el1 = np = 3 19 140 620
el2 = np = 3 154 3 879
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1081 and 20
el1 = np = 3 19 140 620
el2 = np = 3 3 31 165
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1082 and 1725
el1 = np = 3 141 18 621
el2 = np = 3 152 153 908
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1082 and 1760
el1 = np = 3 141 18 621
el2 = np = 3 154 923 153
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1083 and 21
el1 = np = 3 143 144 622
el2 = np = 3 31 32 166
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1083 and 23
el1 = np = 3 143 144 622
el2 = np = 3 32 33 167
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1084 and 21
el1 = np = 3 144 19 623
el2 = np = 3 31 32 166
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1084 and 20
el1 = np = 3 144 19 623
el2 = np = 3 3 31 165
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1085 and 21
el1 = np = 3 144 623 622
el2 = np = 3 31 32 166
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1086 and 1721
el1 = np = 3 17 99 603
el2 = np = 3 150 151 906
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1087 and 25
el1 = np = 3 75 8 624
el2 = np = 3 33 34 168
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1097 and 20
el1 = np = 3 19 620 623
el2 = np = 3 3 31 165
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1103 and 1725
el1 = np = 3 18 17 603
el2 = np = 3 152 153 908
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1103 and 1721
el1 = np = 3 18 17 603
el2 = np = 3 150 151 906
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1104 and 23
el1 = np = 3 142 143 630
el2 = np = 3 32 33 167
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1104 and 25
el1 = np = 3 142 143 630
el2 = np = 3 33 34 168
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1107 and 23
el1 = np = 3 630 143 622
el2 = np = 3 32 33 167
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1177 and 1676
el1 = np = 3 620 140 670
el2 = np = 3 154 3 879
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1179 and 1760
el1 = np = 3 141 621 672
el2 = np = 3 154 923 153
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1180 and 1760
el1 = np = 3 140 141 672
el2 = np = 3 154 923 153
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1180 and 1676
el1 = np = 3 140 141 672
el2 = np = 3 154 3 879
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1181 and 25
el1 = np = 3 8 142 624
el2 = np = 3 33 34 168
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1185 and 1725
el1 = np = 3 18 603 675
el2 = np = 3 152 153 908
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1200 and 1725
el1 = np = 3 18 675 621
el2 = np = 3 152 153 908
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1332 and 1676
el1 = np = 3 670 140 672
el2 = np = 3 154 3 879
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1334 and 25
el1 = np = 3 624 142 630
el2 = np = 3 33 34 168
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1676 and 1332
el1 = np = 3 154 3 879
el2 = np = 3 670 140 672
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1676 and 1180
el1 = np = 3 154 3 879
el2 = np = 3 140 141 672
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1676 and 1177
el1 = np = 3 154 3 879
el2 = np = 3 620 140 670
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1676 and 1081
el1 = np = 3 154 3 879
el2 = np = 3 19 140 620
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1684 and 210
el1 = np = 3 40 4 145
el2 = np = 3 77 78 251
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1684 and 316
el1 = np = 3 40 4 145
el2 = np = 3 77 251 226
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1684 and 161
el1 = np = 3 40 4 145
el2 = np = 3 76 77 226
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1711 and 162
el1 = np = 3 145 9 900
el2 = np = 3 76 226 225
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1711 and 161
el1 = np = 3 145 9 900
el2 = np = 3 76 77 226
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1713 and 167
el1 = np = 3 146 147 902
el2 = np = 3 81 10 230
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1713 and 213
el1 = np = 3 146 147 902
el2 = np = 3 229 81 230
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1713 and 166
el1 = np = 3 146 147 902
el2 = np = 3 11 81 229
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1715 and 169
el1 = np = 3 147 148 903
el2 = np = 3 82 11 231
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1715 and 171
el1 = np = 3 147 148 903
el2 = np = 3 11 229 231
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1715 and 166
el1 = np = 3 147 148 903
el2 = np = 3 11 81 229
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1717 and 168
el1 = np = 3 148 149 904
el2 = np = 3 12 82 99
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1717 and 170
el1 = np = 3 148 149 904
el2 = np = 3 99 82 231
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1717 and 169
el1 = np = 3 148 149 904
el2 = np = 3 82 11 231
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1719 and 168
el1 = np = 3 149 150 905
el2 = np = 3 12 82 99
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1719 and 1053
el1 = np = 3 149 150 905
el2 = np = 3 12 99 17
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1721 and 1053
el1 = np = 3 150 151 906
el2 = np = 3 12 99 17
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1721 and 1086
el1 = np = 3 150 151 906
el2 = np = 3 17 99 603
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1721 and 1103
el1 = np = 3 150 151 906
el2 = np = 3 18 17 603
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1725 and 1185
el1 = np = 3 152 153 908
el2 = np = 3 18 603 675
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1725 and 1200
el1 = np = 3 152 153 908
el2 = np = 3 18 675 621
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1725 and 1082
el1 = np = 3 152 153 908
el2 = np = 3 141 18 621
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1725 and 1103
el1 = np = 3 152 153 908
el2 = np = 3 18 17 603
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1760 and 1082
el1 = np = 3 154 923 153
el2 = np = 3 141 18 621
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1760 and 1179
el1 = np = 3 154 923 153
el2 = np = 3 141 621 672
layer1 = 1
layer2 = 1
 WARNING: Intersecting elements 1760 and 1180
el1 = np = 3 154 923 153
el2 = np = 3 140 141 672
layer1 = 1
layer2 = 1
Num. of mesh nodes     generated: 1114
Num. of mesh triangles generated: 2190
Mesh was generated.
Mesh generation succeeded.

which seems to work (I visualized the mesh and it seems more or less alright).

joachim · February 20, 2024, 12:44pm

can you visualize the whole geometry, and only the critical face in the classical netgen gui, or in webgui ?
This should already trigger the triangulation.

nicoloc · February 20, 2024, 1:03pm

I can see the whole geometry in netgen but probably the critical face is missing (it’s a very small face so I’m not able to tell if it’s there or not with certainty).
These are the logs:

netgen
NETGEN-6.2.2401
Developed by Joachim Schoeberl at
2010-xxxx Vienna University of Technology
2006-2010 RWTH Aachen University
1996-2006 Johannes Kepler University Linz
Including OpenCascade geometry kernel
optfile ./ng.opt does not exist - using default values
togl-version : 2
OCC module loaded
loading ngsolve library
NGSolve-6.2.2401
Using Lapack
Including sparse direct solver UMFPACK
Running parallel using 24 thread(s)
 Load STEP geometry file /tmp/my_part.step
Preparing visualization (deflection = 0.01) ... done
cannot visualize face 4

joachim · February 22, 2024, 5:51am

Hi Nicoloc,

can you try the latest update ?
We are calling the occ-triangulation with IMeshTools parameters now, which seems to be more robust. We also check the error-code, now.

github.com

NGSolve/netgen/blob/4ff7a2261b6a4b3eaf20d86e636427e5afa7afb0/libsrc/occ/occgenmesh.cpp#L672


      
          TopLoc_Location loc;

          Handle(Geom_Surface) surf = BRep_Tool::Surface (face);

          Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);

          

          if (triangulation.IsNull())

            {

              BRepTools::Clean (geom.shape);

              // BRepMesh_IncrementalMesh (geom.shape, 0.01, true);

          

              // https://dev.opencascade.org/doc/overview/html/occt_user_guides__mesh.html

              IMeshTools_Parameters aMeshParams;

              aMeshParams.Deflection               = 0.01;

              aMeshParams.Angle                    = 0.5;

              aMeshParams.Relative                 = Standard_False;

              aMeshParams.InParallel               = Standard_True;

              aMeshParams.MinSize                  = Precision::Confusion();

              aMeshParams.InternalVerticesMode     = Standard_True;

              aMeshParams.ControlSurfaceDeflection = Standard_True;

              

              BRepMesh_IncrementalMesh aMesher (geom.shape, aMeshParams);

              const Standard_Integer aStatus = aMesher.GetStatusFlags();

Joachim

nicoloc · February 22, 2024, 9:16am

That’s very nice, thanks for the change! I tested it on 2 problematic step files. One works now, the other one still gives some issues, I’ll investigate it better!

If I can be more picky, it would be great to know what is the problematic face in the exception, i.e., something like:

std::string errorMessage = "OCC-Triangulation could not be built for face " + std::to_string(i);
throw Exception(errorMessage);

joachim · February 22, 2024, 9:44am

super, good to here it goes into the right direction !

Does the second example throw an exception from the same part of the code ?
The modification is possibly needed at other places as well.

We can only report what we get from OCCT, I only see the status code, but not the face number.
The function OCCSetLocalMeshSize is certainly not the first place for geometry debugging.

Instead, you can loop over the faces in Python, and call face.Triangulation() for each face,
or better the attribute face._webgui_data.

Joachim

nicoloc · February 22, 2024, 9:56am

Does the second example throw an exception from the same part of the code ?

Unfortunately, yes. However, if I remember correctly it used to work with my hacky modification in the past days. I will spend some time trying to understand what changed.

We can only report what we get from OCCT, I only see the status code, but not the face number.
The function OCCSetLocalMeshSize is certainly not the first place for geometry debugging.

I thought that the i variable in the outer for loop would refer to the face number, but I was probably wrong?

github.com

NGSolve/netgen/blob/4ff7a2261b6a4b3eaf20d86e636427e5afa7afb0/libsrc/occ/occgenmesh.cpp#L657


      
                  mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature), mparam), layer);

                }

            }

          

          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));

              int layer = OCCGeometry::GetProperties(face).layer;

              TopLoc_Location loc;

              Handle(Geom_Surface) surf = BRep_Tool::Surface (face);

              Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);

          

              if (triangulation.IsNull())

                {

Curiously, when I checked it out in python the i in the for loop where the exception occurred was exactly referring to the problematic face

Thanks again for your help!