1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include "GmshConfig.h"
7 #include "GmshMessage.h"
8 #include "GModelIO_OCC.h"
9 #include "Context.h"
10 #include "OCCVertex.h"
11 #include "OCCEdge.h"
12 #include "OCCFace.h"
13 #include "OCCRegion.h"
14 #include "MElement.h"
15 #include "MLine.h"
16 #include "OpenFile.h"
17 #include "StringUtils.h"
18 #include "ExtrudeParams.h"
19
20 #if defined(HAVE_OCC)
21
22 #include <BRepAlgoAPI_Common.hxx>
23 #include <BRepAlgoAPI_Cut.hxx>
24 #include <BRepAlgoAPI_Fuse.hxx>
25 #include <BRepAlgoAPI_Section.hxx>
26 #include <BRepBndLib.hxx>
27 #include <BRepBuilderAPI_Copy.hxx>
28 #include <BRepBuilderAPI_GTransform.hxx>
29 #include <BRepBuilderAPI_MakeEdge.hxx>
30 #include <BRepBuilderAPI_MakeFace.hxx>
31 #include <BRepBuilderAPI_MakeShell.hxx>
32 #include <BRepBuilderAPI_MakeSolid.hxx>
33 #include <BRepBuilderAPI_MakeVertex.hxx>
34 #include <BRepBuilderAPI_MakeWire.hxx>
35 #include <BRepBuilderAPI_Transform.hxx>
36 #include <BRepCheck_Analyzer.hxx>
37 #include <BRepExtrema_DistShapeShape.hxx>
38 #include <BRepFill_CurveConstraint.hxx>
39 #include <BRepFilletAPI_MakeChamfer.hxx>
40 #include <BRepFilletAPI_MakeFillet.hxx>
41 #include <BRepGProp.hxx>
42 #include <BRepLib.hxx>
43 #include <BRepOffsetAPI_MakeFilling.hxx>
44 #include <BRepOffsetAPI_MakePipe.hxx>
45 #include <BRepOffsetAPI_MakeThickSolid.hxx>
46 #include <BRepOffsetAPI_Sewing.hxx>
47 #include <BRepOffsetAPI_ThruSections.hxx>
48 #include <BRepPrimAPI_MakeBox.hxx>
49 #include <BRepPrimAPI_MakeCone.hxx>
50 #include <BRepPrimAPI_MakeCylinder.hxx>
51 #include <BRepPrimAPI_MakePrism.hxx>
52 #include <BRepPrimAPI_MakeRevol.hxx>
53 #include <BRepPrimAPI_MakeSphere.hxx>
54 #include <BRepPrimAPI_MakeTorus.hxx>
55 #include <BRepPrimAPI_MakeWedge.hxx>
56 #include <BRepTools.hxx>
57 #include <BRep_Tool.hxx>
58 #include <Bnd_Box.hxx>
59 #include <ElCLib.hxx>
60 #include <GProp_GProps.hxx>
61 #include <Geom2d_Curve.hxx>
62 #include <GeomAPI_Interpolate.hxx>
63 #include <GeomFill_BSplineCurves.hxx>
64 #include <GeomFill_BezierCurves.hxx>
65 #include <GeomProjLib.hxx>
66 #include <Geom_BSplineCurve.hxx>
67 #include <Geom_BSplineSurface.hxx>
68 #include <Geom_BezierCurve.hxx>
69 #include <Geom_BezierSurface.hxx>
70 #include <Geom_Circle.hxx>
71 #include <Geom_Ellipse.hxx>
72 #include <Geom_Plane.hxx>
73 #include <Geom_Surface.hxx>
74 #include <Geom_TrimmedCurve.hxx>
75 #include <IGESControl_Reader.hxx>
76 #include <IGESControl_Writer.hxx>
77 #include <Interface_Static.hxx>
78 #include <Poly_PolygonOnTriangulation.hxx>
79 #include <Poly_Triangle.hxx>
80 #include <Poly_Triangulation.hxx>
81 #include <ProjLib_ProjectedCurve.hxx>
82 #include <STEPControl_Reader.hxx>
83 #include <STEPControl_Writer.hxx>
84 #include <ShapeBuild_ReShape.hxx>
85 #include <ShapeFix_FixSmallFace.hxx>
86 #include <ShapeFix_Shape.hxx>
87 #include <ShapeFix_Wireframe.hxx>
88 #include <Standard_Version.hxx>
89 #include <TColStd_Array1OfInteger.hxx>
90 #include <TColStd_Array1OfReal.hxx>
91 #include <TColStd_Array2OfReal.hxx>
92 #include <TColgp_Array1OfPnt.hxx>
93 #include <TColgp_Array1OfPnt2d.hxx>
94 #include <TColgp_Array2OfPnt.hxx>
95 #include <TColgp_HArray1OfPnt.hxx>
96 #include <TopExp.hxx>
97 #include <TopExp_Explorer.hxx>
98 #include <TopTools_DataMapIteratorOfDataMapOfIntegerShape.hxx>
99 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
100 #include <TopTools_ListIteratorOfListOfShape.hxx>
101 #include <TopoDS.hxx>
102 #include <gce_MakeCirc.hxx>
103 #include <gce_MakeElips.hxx>
104 #include <gce_MakePln.hxx>
105 #include <numeric>
106 #include <utility>
107
108 #include "OCCAttributes.h"
109
110 #if OCC_VERSION_HEX < 0x060900
111 #error "Gmsh requires OpenCASCADE >= 6.9"
112 #endif
113
114 #if OCC_VERSION_HEX > 0x070300
115 #include <BRepMesh_IncrementalMesh.hxx>
116 #else
117 #include <BRepMesh_FastDiscret.hxx>
118 #endif
119
120 #if OCC_VERSION_HEX < 0x070400
121 #include <ShapeUpgrade_UnifySameDomain.hxx>
122 #endif
123
124 #if defined(HAVE_OCC_CAF)
125 #include <IGESCAFControl_Reader.hxx>
126 #include <Quantity_Color.hxx>
127 #include <STEPCAFControl_Reader.hxx>
128 #include <TDF_ChildIterator.hxx>
129 #include <TDF_Tool.hxx>
130 #include <TDataStd_Name.hxx>
131 #include <TDocStd_Document.hxx>
132 #include <XCAFApp_Application.hxx>
133 #include <XCAFDoc_Color.hxx>
134 #include <XCAFDoc_ColorTool.hxx>
135 #include <XCAFDoc_DocumentTool.hxx>
136 #include <XCAFDoc_Location.hxx>
137 #include <XCAFDoc_MaterialTool.hxx>
138 #include <XCAFDoc_ShapeTool.hxx>
139 #endif
140
OCC_Internals()141 OCC_Internals::OCC_Internals()
142 {
143 for(int i = 0; i < 6; i++) _maxTag[i] = 0;
144 _changed = true;
145 _attributes = new OCCAttributesRTree(CTX::instance()->geom.tolerance);
146 }
147
~OCC_Internals()148 OCC_Internals::~OCC_Internals() { delete _attributes; }
149
reset()150 void OCC_Internals::reset()
151 {
152 _attributes->clear();
153 _somap.Clear();
154 _shmap.Clear();
155 _fmap.Clear();
156 _wmap.Clear();
157 _emap.Clear();
158 _vmap.Clear();
159 _unbind();
160 }
161
setMaxTag(int dim,int val)162 void OCC_Internals::setMaxTag(int dim, int val)
163 {
164 if(dim < -2 || dim > 3) return;
165 _maxTag[dim + 2] = std::max(_maxTag[dim + 2], val);
166 }
167
getMaxTag(int dim) const168 int OCC_Internals::getMaxTag(int dim) const
169 {
170 if(dim < -2 || dim > 3) return 0;
171 return _maxTag[dim + 2];
172 }
173
_recomputeMaxTag(int dim)174 void OCC_Internals::_recomputeMaxTag(int dim)
175 {
176 if(dim < -2 || dim > 3) return;
177 _maxTag[dim + 2] = 0;
178 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp;
179 switch(dim) {
180 case 0: exp.Initialize(_tagVertex); break;
181 case 1: exp.Initialize(_tagEdge); break;
182 case 2: exp.Initialize(_tagFace); break;
183 case 3: exp.Initialize(_tagSolid); break;
184 case -1: exp.Initialize(_tagWire); break;
185 case -2: exp.Initialize(_tagShell); break;
186 }
187 for(; exp.More(); exp.Next())
188 _maxTag[dim + 2] = std::max(_maxTag[dim + 2], exp.Key());
189 }
190
_bind(const TopoDS_Vertex & vertex,int tag,bool recursive)191 void OCC_Internals::_bind(const TopoDS_Vertex &vertex, int tag, bool recursive)
192 {
193 if(vertex.IsNull()) return;
194 if(_vertexTag.IsBound(vertex)) {
195 if(_vertexTag.Find(vertex) != tag) {
196 Msg::Info("Cannot bind existing OpenCASCADE point %d to second tag %d",
197 _vertexTag.Find(vertex), tag);
198 }
199 }
200 else {
201 if(_tagVertex.IsBound(tag)) {
202 // this leaves the old vertex bound in _vertexTag, but we cannot remove it
203 Msg::Info("Rebinding OpenCASCADE point %d", tag);
204 }
205 _vertexTag.Bind(vertex, tag);
206 _tagVertex.Bind(tag, vertex);
207 setMaxTag(0, tag);
208 _changed = true;
209 _attributes->insert(new OCCAttributes(0, vertex));
210 }
211 }
212
_bind(const TopoDS_Edge & edge,int tag,bool recursive)213 void OCC_Internals::_bind(const TopoDS_Edge &edge, int tag, bool recursive)
214 {
215 if(edge.IsNull()) return;
216 if(_edgeTag.IsBound(edge)) {
217 if(_edgeTag.Find(edge) != tag) {
218 Msg::Info("Cannot bind existing OpenCASCADE curve %d to second tag %d",
219 _edgeTag.Find(edge), tag);
220 }
221 }
222 else {
223 if(_tagEdge.IsBound(tag)) {
224 // this leaves the old edge bound in _edgeTag, but we cannot remove it
225 Msg::Info("Rebinding OpenCASCADE curve %d", tag);
226 }
227 _edgeTag.Bind(edge, tag);
228 _tagEdge.Bind(tag, edge);
229 setMaxTag(1, tag);
230 _changed = true;
231 _attributes->insert(new OCCAttributes(1, edge));
232 }
233 if(recursive) {
234 TopExp_Explorer exp0;
235 for(exp0.Init(edge, TopAbs_VERTEX); exp0.More(); exp0.Next()) {
236 TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
237 if(!_vertexTag.IsBound(vertex)) {
238 int t = getMaxTag(0) + 1;
239 _bind(vertex, t, recursive);
240 }
241 }
242 }
243 }
244
_bind(const TopoDS_Wire & wire,int tag,bool recursive)245 void OCC_Internals::_bind(const TopoDS_Wire &wire, int tag, bool recursive)
246 {
247 if(wire.IsNull()) return;
248 if(_wireTag.IsBound(wire)) {
249 if(_wireTag.Find(wire) != tag) {
250 Msg::Info("Cannot bind existing OpenCASCADE wire %d to second tag %d",
251 _wireTag.Find(wire), tag);
252 }
253 }
254 else {
255 if(_tagWire.IsBound(tag)) {
256 // this leaves the old wire bound in _wireTag, but we cannot remove it
257 Msg::Info("Rebinding OpenCASCADE wire %d", tag);
258 }
259 _wireTag.Bind(wire, tag);
260 _tagWire.Bind(tag, wire);
261 setMaxTag(-1, tag);
262 _changed = true;
263 }
264 if(recursive) {
265 TopExp_Explorer exp0;
266 for(exp0.Init(wire, TopAbs_EDGE); exp0.More(); exp0.Next()) {
267 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
268 if(!_edgeTag.IsBound(edge)) {
269 int t = getMaxTag(1) + 1;
270 _bind(edge, t, recursive);
271 }
272 }
273 }
274 }
275
_bind(const TopoDS_Face & face,int tag,bool recursive)276 void OCC_Internals::_bind(const TopoDS_Face &face, int tag, bool recursive)
277 {
278 if(face.IsNull()) return;
279 if(_faceTag.IsBound(face)) {
280 if(_faceTag.Find(face) != tag) {
281 Msg::Info("Cannot bind existing OpenCASCADE surface %d to second tag %d",
282 _faceTag.Find(face), tag);
283 }
284 }
285 else {
286 if(_tagFace.IsBound(tag)) {
287 // this leaves the old face bound in _faceTag, but we cannot remove it
288 Msg::Info("Rebinding OpenCASCADE surface %d", tag);
289 }
290 _faceTag.Bind(face, tag);
291 _tagFace.Bind(tag, face);
292 setMaxTag(2, tag);
293 _changed = true;
294 _attributes->insert(new OCCAttributes(2, face));
295 }
296 if(recursive) {
297 TopExp_Explorer exp0;
298 for(exp0.Init(face, TopAbs_WIRE); exp0.More(); exp0.Next()) {
299 TopoDS_Wire wire = TopoDS::Wire(exp0.Current());
300 if(!_wireTag.IsBound(wire)) {
301 int t = getMaxTag(-1) + 1;
302 _bind(wire, t, recursive);
303 }
304 }
305 for(exp0.Init(face, TopAbs_EDGE); exp0.More(); exp0.Next()) {
306 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
307 if(!_edgeTag.IsBound(edge)) {
308 int t = getMaxTag(1) + 1;
309 _bind(edge, t, recursive);
310 }
311 }
312 }
313 }
314
_bind(const TopoDS_Shell & shell,int tag,bool recursive)315 void OCC_Internals::_bind(const TopoDS_Shell &shell, int tag, bool recursive)
316 {
317 if(shell.IsNull()) return;
318 if(_shellTag.IsBound(shell)) {
319 if(_shellTag.Find(shell) != tag) {
320 Msg::Info("Cannot bind existing OpenCASCADE shell %d to second tag %d",
321 _shellTag.Find(shell), tag);
322 }
323 }
324 else {
325 if(_tagShell.IsBound(tag)) {
326 // this leaves the old shell bound in _faceTag, but we cannot remove it
327 Msg::Info("Rebinding OpenCASCADE shell %d", tag);
328 }
329 _shellTag.Bind(shell, tag);
330 _tagShell.Bind(tag, shell);
331 setMaxTag(-2, tag);
332 _changed = true;
333 }
334 if(recursive) {
335 TopExp_Explorer exp0;
336 for(exp0.Init(shell, TopAbs_FACE); exp0.More(); exp0.Next()) {
337 TopoDS_Face face = TopoDS::Face(exp0.Current());
338 if(!_faceTag.IsBound(face)) {
339 int t = getMaxTag(2) + 1;
340 _bind(face, t, recursive);
341 }
342 }
343 }
344 }
345
_bind(const TopoDS_Solid & solid,int tag,bool recursive)346 void OCC_Internals::_bind(const TopoDS_Solid &solid, int tag, bool recursive)
347 {
348 if(solid.IsNull()) return;
349 if(_solidTag.IsBound(solid)) {
350 if(_solidTag.Find(solid) != tag) {
351 Msg::Info("Cannot bind existing OpenCASCADE volume %d to second tag %d",
352 _solidTag.Find(solid), tag);
353 }
354 }
355 else {
356 if(_tagSolid.IsBound(tag)) {
357 // this leaves the old solid bound in _faceTag, but we cannot remove it
358 Msg::Info("Rebinding OpenCASCADE volume %d", tag);
359 }
360 _solidTag.Bind(solid, tag);
361 _tagSolid.Bind(tag, solid);
362 setMaxTag(3, tag);
363 _changed = true;
364 _attributes->insert(new OCCAttributes(3, solid));
365 }
366 if(recursive) {
367 TopExp_Explorer exp0;
368 for(exp0.Init(solid, TopAbs_SHELL); exp0.More(); exp0.Next()) {
369 TopoDS_Shell shell = TopoDS::Shell(exp0.Current());
370 if(!_shellTag.IsBound(shell)) {
371 int t = getMaxTag(-2) + 1;
372 _bind(shell, t, recursive);
373 }
374 }
375 for(exp0.Init(solid, TopAbs_FACE); exp0.More(); exp0.Next()) {
376 TopoDS_Face face = TopoDS::Face(exp0.Current());
377 if(!_faceTag.IsBound(face)) {
378 int t = getMaxTag(2) + 1;
379 _bind(face, t, recursive);
380 }
381 }
382 }
383 }
384
_bind(TopoDS_Shape shape,int dim,int tag,bool recursive)385 void OCC_Internals::_bind(TopoDS_Shape shape, int dim, int tag, bool recursive)
386 {
387 switch(dim) {
388 case 0: _bind(TopoDS::Vertex(shape), tag, recursive); break;
389 case 1: _bind(TopoDS::Edge(shape), tag, recursive); break;
390 case 2: _bind(TopoDS::Face(shape), tag, recursive); break;
391 case 3: _bind(TopoDS::Solid(shape), tag, recursive); break;
392 case -1: _bind(TopoDS::Wire(shape), tag, recursive); break;
393 case -2: _bind(TopoDS::Shell(shape), tag, recursive); break;
394 default: break;
395 }
396 }
397
_unbind(const TopoDS_Vertex & vertex,int tag,bool recursive)398 void OCC_Internals::_unbind(const TopoDS_Vertex &vertex, int tag, bool recursive)
399 {
400 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp0(_tagEdge);
401 for(; exp0.More(); exp0.Next()) {
402 TopoDS_Edge edge = TopoDS::Edge(exp0.Value());
403 TopExp_Explorer exp1;
404 for(exp1.Init(edge, TopAbs_VERTEX); exp1.More(); exp1.Next()) {
405 if(exp1.Current().IsSame(vertex)) return;
406 }
407 }
408 std::pair<int, int> dimTag(0, tag);
409 if(_toPreserve.find(dimTag) != _toPreserve.end()) return;
410 _vertexTag.UnBind(vertex);
411 _tagVertex.UnBind(tag);
412 _toRemove.insert(dimTag);
413 _recomputeMaxTag(0);
414 _changed = true;
415 }
416
_unbind(const TopoDS_Edge & edge,int tag,bool recursive)417 void OCC_Internals::_unbind(const TopoDS_Edge &edge, int tag, bool recursive)
418 {
419 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp2(_tagFace);
420 for(; exp2.More(); exp2.Next()) {
421 TopoDS_Face face = TopoDS::Face(exp2.Value());
422 TopExp_Explorer exp1;
423 for(exp1.Init(face, TopAbs_EDGE); exp1.More(); exp1.Next()) {
424 if(exp1.Current().IsSame(edge)) return;
425 }
426 }
427 std::pair<int, int> dimTag(1, tag);
428 if(_toPreserve.find(dimTag) != _toPreserve.end()) return;
429 _edgeTag.UnBind(edge);
430 _tagEdge.UnBind(tag);
431 _toRemove.insert(dimTag);
432 _recomputeMaxTag(1);
433 if(recursive) {
434 TopExp_Explorer exp0;
435 for(exp0.Init(edge, TopAbs_VERTEX); exp0.More(); exp0.Next()) {
436 TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
437 if(_vertexTag.IsBound(vertex)) {
438 int t = _vertexTag.Find(vertex);
439 _unbind(vertex, t, recursive);
440 }
441 }
442 }
443 _changed = true;
444 }
445
_unbind(const TopoDS_Wire & wire,int tag,bool recursive)446 void OCC_Internals::_unbind(const TopoDS_Wire &wire, int tag, bool recursive)
447 {
448 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp0(_tagFace);
449 for(; exp0.More(); exp0.Next()) {
450 TopoDS_Face face = TopoDS::Face(exp0.Value());
451 TopExp_Explorer exp1;
452 for(exp1.Init(face, TopAbs_WIRE); exp1.More(); exp1.Next()) {
453 if(exp1.Current().IsSame(wire)) return;
454 }
455 }
456 std::pair<int, int> dimTag(-1, tag);
457 if(_toPreserve.find(dimTag) != _toPreserve.end()) return;
458 _wireTag.UnBind(wire);
459 _tagWire.UnBind(tag);
460 _toRemove.insert(dimTag);
461 _recomputeMaxTag(-1);
462 if(recursive) {
463 TopExp_Explorer exp0;
464 for(exp0.Init(wire, TopAbs_EDGE); exp0.More(); exp0.Next()) {
465 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
466 if(_edgeTag.IsBound(edge)) {
467 int t = _edgeTag.Find(edge);
468 _unbind(edge, t, recursive);
469 }
470 }
471 }
472 _changed = true;
473 }
474
_unbind(const TopoDS_Face & face,int tag,bool recursive)475 void OCC_Internals::_unbind(const TopoDS_Face &face, int tag, bool recursive)
476 {
477 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp2(_tagSolid);
478 for(; exp2.More(); exp2.Next()) {
479 TopoDS_Solid solid = TopoDS::Solid(exp2.Value());
480 TopExp_Explorer exp1;
481 for(exp1.Init(solid, TopAbs_FACE); exp1.More(); exp1.Next()) {
482 if(exp1.Current().IsSame(face)) return;
483 }
484 }
485 std::pair<int, int> dimTag(2, tag);
486 if(_toPreserve.find(dimTag) != _toPreserve.end()) return;
487 _faceTag.UnBind(face);
488 _tagFace.UnBind(tag);
489 _toRemove.insert(dimTag);
490 _recomputeMaxTag(2);
491 if(recursive) {
492 TopExp_Explorer exp0;
493 for(exp0.Init(face, TopAbs_WIRE); exp0.More(); exp0.Next()) {
494 TopoDS_Wire wire = TopoDS::Wire(exp0.Current());
495 if(_wireTag.IsBound(wire)) {
496 int t = _wireTag.Find(wire);
497 _unbind(wire, t, recursive);
498 }
499 }
500 for(exp0.Init(face, TopAbs_EDGE); exp0.More(); exp0.Next()) {
501 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
502 if(_edgeTag.IsBound(edge)) {
503 int t = _edgeTag.Find(edge);
504 _unbind(edge, t, recursive);
505 }
506 }
507 }
508 _changed = true;
509 }
510
_unbind(const TopoDS_Shell & shell,int tag,bool recursive)511 void OCC_Internals::_unbind(const TopoDS_Shell &shell, int tag, bool recursive)
512 {
513 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp0(_tagSolid);
514 for(; exp0.More(); exp0.Next()) {
515 TopoDS_Solid solid = TopoDS::Solid(exp0.Value());
516 TopExp_Explorer exp1;
517 for(exp1.Init(solid, TopAbs_SHELL); exp1.More(); exp1.Next()) {
518 if(exp1.Current().IsSame(shell)) return;
519 }
520 }
521 std::pair<int, int> dimTag(-2, tag);
522 if(_toPreserve.find(dimTag) != _toPreserve.end()) return;
523 _shellTag.UnBind(shell);
524 _tagShell.UnBind(tag);
525 _toRemove.insert(dimTag);
526 _recomputeMaxTag(-2);
527 if(recursive) {
528 TopExp_Explorer exp0;
529 for(exp0.Init(shell, TopAbs_FACE); exp0.More(); exp0.Next()) {
530 TopoDS_Face face = TopoDS::Face(exp0.Current());
531 if(_faceTag.IsBound(face)) {
532 int t = _faceTag.Find(face);
533 _unbind(face, t, recursive);
534 }
535 }
536 }
537 _changed = true;
538 }
539
_unbind(const TopoDS_Solid & solid,int tag,bool recursive)540 void OCC_Internals::_unbind(const TopoDS_Solid &solid, int tag, bool recursive)
541 {
542 std::pair<int, int> dimTag(3, tag);
543 if(_toPreserve.find(dimTag) != _toPreserve.end()) return;
544 _solidTag.UnBind(solid);
545 _tagSolid.UnBind(tag);
546 _toRemove.insert(dimTag);
547 _recomputeMaxTag(3);
548 if(recursive) {
549 TopExp_Explorer exp0;
550 for(exp0.Init(solid, TopAbs_SHELL); exp0.More(); exp0.Next()) {
551 TopoDS_Shell shell = TopoDS::Shell(exp0.Current());
552 if(_shellTag.IsBound(shell)) {
553 int t = _shellTag.Find(shell);
554 _unbind(shell, t, recursive);
555 }
556 }
557 for(exp0.Init(solid, TopAbs_FACE); exp0.More(); exp0.Next()) {
558 TopoDS_Face face = TopoDS::Face(exp0.Current());
559 if(_faceTag.IsBound(face)) {
560 int t = _faceTag.Find(face);
561 _unbind(face, t, recursive);
562 }
563 }
564 }
565 _changed = true;
566 }
567
_unbind(TopoDS_Shape shape,int dim,int tag,bool recursive)568 void OCC_Internals::_unbind(TopoDS_Shape shape, int dim, int tag, bool recursive)
569 {
570 switch(dim) {
571 case 0: _unbind(TopoDS::Vertex(shape), tag, recursive); break;
572 case 1: _unbind(TopoDS::Edge(shape), tag, recursive); break;
573 case 2: _unbind(TopoDS::Face(shape), tag, recursive); break;
574 case 3: _unbind(TopoDS::Solid(shape), tag, recursive); break;
575 case -1: _unbind(TopoDS::Wire(shape), tag, recursive); break;
576 case -2: _unbind(TopoDS::Shell(shape), tag, recursive); break;
577 default: break;
578 }
579 }
580
_unbindWithoutChecks(TopoDS_Shape shape)581 void OCC_Internals::_unbindWithoutChecks(TopoDS_Shape shape)
582 {
583 TopExp_Explorer exp0;
584 for(exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()) {
585 TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
586 if(_solidTag.IsBound(solid)) {
587 int tag = _solidTag.Find(solid);
588 _solidTag.UnBind(solid);
589 _tagSolid.UnBind(tag);
590 _toRemove.insert(std::make_pair(3, tag));
591 _changed = true;
592 }
593 }
594 for(exp0.Init(shape, TopAbs_SHELL); exp0.More(); exp0.Next()) {
595 TopoDS_Shell shell = TopoDS::Shell(exp0.Current());
596 if(_shellTag.IsBound(shell)) {
597 int tag = _shellTag.Find(shell);
598 _shellTag.UnBind(shell);
599 _tagShell.UnBind(tag);
600 _toRemove.insert(std::make_pair(-2, tag));
601 _changed = true;
602 }
603 }
604 for(exp0.Init(shape, TopAbs_FACE); exp0.More(); exp0.Next()) {
605 TopoDS_Face face = TopoDS::Face(exp0.Current());
606 if(_faceTag.IsBound(face)) {
607 int tag = _faceTag.Find(face);
608 _faceTag.UnBind(face);
609 _tagFace.UnBind(tag);
610 _toRemove.insert(std::make_pair(2, tag));
611 _changed = true;
612 }
613 }
614 for(exp0.Init(shape, TopAbs_WIRE); exp0.More(); exp0.Next()) {
615 TopoDS_Wire wire = TopoDS::Wire(exp0.Current());
616 if(_wireTag.IsBound(wire)) {
617 int tag = _wireTag.Find(wire);
618 _wireTag.UnBind(wire);
619 _tagWire.UnBind(tag);
620 _toRemove.insert(std::make_pair(-1, tag));
621 _changed = true;
622 }
623 }
624 for(exp0.Init(shape, TopAbs_EDGE); exp0.More(); exp0.Next()) {
625 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
626 if(_edgeTag.IsBound(edge)) {
627 int tag = _edgeTag.Find(edge);
628 _edgeTag.UnBind(edge);
629 _tagEdge.UnBind(tag);
630 _toRemove.insert(std::make_pair(1, tag));
631 _changed = true;
632 }
633 }
634 for(exp0.Init(shape, TopAbs_VERTEX); exp0.More(); exp0.Next()) {
635 TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
636 if(_vertexTag.IsBound(vertex)) {
637 int tag = _vertexTag.Find(vertex);
638 _vertexTag.UnBind(vertex);
639 _tagVertex.UnBind(tag);
640 _toRemove.insert(std::make_pair(0, tag));
641 _changed = true;
642 }
643 }
644 }
645
_unbind()646 void OCC_Internals::_unbind()
647 {
648 for(int i = 0; i < 6; i++) _maxTag[i] = 0;
649
650 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp;
651 exp.Initialize(_tagVertex);
652 for(; exp.More(); exp.Next())
653 _toRemove.insert(std::make_pair(0, exp.Key()));
654 exp.Initialize(_tagEdge);
655 for(; exp.More(); exp.Next())
656 _toRemove.insert(std::make_pair(1, exp.Key()));
657 exp.Initialize(_tagFace);
658 for(; exp.More(); exp.Next())
659 _toRemove.insert(std::make_pair(2, exp.Key()));
660 exp.Initialize(_tagSolid);
661 for(; exp.More(); exp.Next())
662 _toRemove.insert(std::make_pair(3, exp.Key()));
663 exp.Initialize(_tagWire);
664 for(; exp.More(); exp.Next())
665 _toRemove.insert(std::make_pair(-1, exp.Key()));
666 exp.Initialize(_tagShell);
667 for(; exp.More(); exp.Next())
668 _toRemove.insert(std::make_pair(-2, exp.Key()));
669
670 _tagVertex.Clear();
671 _tagEdge.Clear();
672 _tagFace.Clear();
673 _tagSolid.Clear();
674 _tagWire.Clear();
675 _tagShell.Clear();
676
677 _vertexTag.Clear();
678 _edgeTag.Clear();
679 _faceTag.Clear();
680 _solidTag.Clear();
681 _wireTag.Clear();
682 _shellTag.Clear();
683
684 _changed = true;
685 }
686
_multiBind(const TopoDS_Shape & shape,int tag,std::vector<std::pair<int,int>> & outDimTags,bool highestDimOnly,bool recursive,bool returnNewOnly)687 void OCC_Internals::_multiBind(const TopoDS_Shape &shape, int tag,
688 std::vector<std::pair<int, int> > &outDimTags,
689 bool highestDimOnly, bool recursive,
690 bool returnNewOnly)
691 {
692 TopExp_Explorer exp0;
693 int count = 0;
694 for(exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()) {
695 TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
696 bool exists = false;
697 int t = tag;
698 if(t <= 0) {
699 if(_solidTag.IsBound(solid)) {
700 t = _solidTag.Find(solid);
701 exists = true;
702 }
703 else
704 t = getMaxTag(3) + 1;
705 }
706 else if(count) {
707 Msg::Error("Cannot bind multiple volumes to single tag %d", t);
708 return;
709 }
710 if(!exists) _bind(solid, t, recursive);
711 if(!exists || !returnNewOnly)
712 outDimTags.push_back(std::make_pair(3, t));
713 count++;
714 }
715 if(highestDimOnly && count) return;
716 for(exp0.Init(shape, TopAbs_FACE); exp0.More(); exp0.Next()) {
717 TopoDS_Face face = TopoDS::Face(exp0.Current());
718 bool exists = false;
719 int t = tag;
720 if(t <= 0) {
721 if(_faceTag.IsBound(face)) {
722 t = _faceTag.Find(face);
723 exists = true;
724 }
725 else
726 t = getMaxTag(2) + 1;
727 }
728 else if(count) {
729 Msg::Error("Cannot bind multiple surfaces to single tag %d", t);
730 return;
731 }
732 if(!exists) _bind(face, t, recursive);
733 if(!exists || !returnNewOnly)
734 outDimTags.push_back(std::make_pair(2, t));
735 count++;
736 }
737 if(highestDimOnly && count) return;
738 for(exp0.Init(shape, TopAbs_EDGE); exp0.More(); exp0.Next()) {
739 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
740 bool exists = false;
741 int t = tag;
742 if(t <= 0) {
743 if(_edgeTag.IsBound(edge)) {
744 t = _edgeTag.Find(edge);
745 exists = true;
746 }
747 else
748 t = getMaxTag(1) + 1;
749 }
750 else if(count) {
751 Msg::Error("Cannot bind multiple curves to single tag %d", t);
752 return;
753 }
754 if(!exists) _bind(edge, t, recursive);
755 if(!exists || !returnNewOnly)
756 outDimTags.push_back(std::make_pair(1, t));
757 count++;
758 }
759 if(highestDimOnly && count) return;
760 for(exp0.Init(shape, TopAbs_VERTEX); exp0.More(); exp0.Next()) {
761 TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
762 bool exists = false;
763 int t = tag;
764 if(t <= 0) {
765 if(_vertexTag.IsBound(vertex)) {
766 t = _vertexTag.Find(vertex);
767 exists = true;
768 }
769 else
770 t = getMaxTag(0) + 1;
771 }
772 else if(count) {
773 Msg::Error("Cannot bind multiple points to single tag %d", t);
774 return;
775 }
776 if(!exists) _bind(vertex, t, recursive);
777 if(!exists || !returnNewOnly)
778 outDimTags.push_back(std::make_pair(0, t));
779 count++;
780 }
781 }
782
_isBound(int dim,int tag)783 bool OCC_Internals::_isBound(int dim, int tag)
784 {
785 switch(dim) {
786 case 0: return _tagVertex.IsBound(tag);
787 case 1: return _tagEdge.IsBound(tag);
788 case 2: return _tagFace.IsBound(tag);
789 case 3: return _tagSolid.IsBound(tag);
790 case -1: return _tagWire.IsBound(tag);
791 case -2: return _tagShell.IsBound(tag);
792 default: return false;
793 }
794 }
795
_isBound(int dim,const TopoDS_Shape & shape)796 bool OCC_Internals::_isBound(int dim, const TopoDS_Shape &shape)
797 {
798 switch(dim) {
799 case 0: return _vertexTag.IsBound(shape);
800 case 1: return _edgeTag.IsBound(shape);
801 case 2: return _faceTag.IsBound(shape);
802 case 3: return _solidTag.IsBound(shape);
803 case -1: return _wireTag.IsBound(shape);
804 case -2: return _shellTag.IsBound(shape);
805 default: return false;
806 }
807 }
808
_find(int dim,int tag)809 TopoDS_Shape OCC_Internals::_find(int dim, int tag)
810 {
811 switch(dim) {
812 case 0: return _tagVertex.Find(tag);
813 case 1: return _tagEdge.Find(tag);
814 case 2: return _tagFace.Find(tag);
815 case 3: return _tagSolid.Find(tag);
816 case -1: return _tagWire.Find(tag);
817 case -2: return _tagShell.Find(tag);
818 default: return TopoDS_Shape();
819 }
820 }
821
_find(int dim,const TopoDS_Shape & shape)822 int OCC_Internals::_find(int dim, const TopoDS_Shape &shape)
823 {
824 switch(dim) {
825 case 0: return _vertexTag.Find(shape);
826 case 1: return _edgeTag.Find(shape);
827 case 2: return _faceTag.Find(shape);
828 case 3: return _solidTag.Find(shape);
829 case -1: return _wireTag.Find(shape);
830 case -2: return _shellTag.Find(shape);
831 default: return -1;
832 }
833 }
834
addVertex(int & tag,double x,double y,double z,double meshSize)835 bool OCC_Internals::addVertex(int &tag, double x, double y, double z,
836 double meshSize)
837 {
838 if(tag >= 0 && _tagVertex.IsBound(tag)) {
839 Msg::Error("OpenCASCADE point with tag %d already exists", tag);
840 return false;
841 }
842 TopoDS_Vertex result;
843 try {
844 gp_Pnt aPnt(x, y, z);
845 BRepBuilderAPI_MakeVertex v(aPnt);
846 v.Build();
847 if(!v.IsDone()) {
848 Msg::Error("Could not create point");
849 return false;
850 }
851 result = v.Vertex();
852 } catch(Standard_Failure &err) {
853 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
854 return false;
855 }
856 if(meshSize > 0 && meshSize < MAX_LC)
857 _attributes->insert(new OCCAttributes(0, result, meshSize));
858 if(tag < 0) tag = getMaxTag(0) + 1;
859 _bind(result, tag, true);
860 return true;
861 }
862
addLine(int & tag,int startTag,int endTag)863 bool OCC_Internals::addLine(int &tag, int startTag, int endTag)
864 {
865 if(tag >= 0 && _tagEdge.IsBound(tag)) {
866 Msg::Error("OpenCASCADE curve with tag %d already exists", tag);
867 return false;
868 }
869 if(!_tagVertex.IsBound(startTag)) {
870 Msg::Error("Unknown OpenCASCADE point with tag %d", startTag);
871 return false;
872 }
873 if(!_tagVertex.IsBound(endTag)) {
874 Msg::Error("Unknown OpenCASCADE point with tag %d", endTag);
875 return false;
876 }
877 if(startTag == endTag) {
878 Msg::Error("Start and end points of line should be different");
879 return false;
880 }
881 TopoDS_Edge result;
882 try {
883 TopoDS_Vertex start = TopoDS::Vertex(_tagVertex.Find(startTag));
884 TopoDS_Vertex end = TopoDS::Vertex(_tagVertex.Find(endTag));
885 BRepBuilderAPI_MakeEdge e(start, end);
886 e.Build();
887 if(!e.IsDone()) {
888 Msg::Error("Could not create line");
889 return false;
890 }
891 result = e.Edge();
892 } catch(Standard_Failure &err) {
893 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
894 return false;
895 }
896 if(tag < 0) tag = getMaxTag(1) + 1;
897 _bind(result, tag, true);
898 return true;
899 }
900
addLine(int & tag,const std::vector<int> & pointTags)901 bool OCC_Internals::addLine(int &tag, const std::vector<int> &pointTags)
902 {
903 if(pointTags.size() == 2) return addLine(tag, pointTags[0], pointTags[1]);
904
905 // FIXME: if tag < 0 we could create multiple lines
906 Msg::Error("OpenCASCADE polyline currently not supported");
907 return false;
908 }
909
addCircleArc(int & tag,int startTag,int centerTag,int endTag)910 bool OCC_Internals::addCircleArc(int &tag, int startTag, int centerTag,
911 int endTag)
912 {
913 if(tag >= 0 && _tagEdge.IsBound(tag)) {
914 Msg::Error("OpenCASCADE curve with tag %d already exists", tag);
915 return false;
916 }
917 if(!_tagVertex.IsBound(startTag)) {
918 Msg::Error("Unknown OpenCASCADE point with tag %d", startTag);
919 return false;
920 }
921 if(!_tagVertex.IsBound(centerTag)) {
922 Msg::Error("Unknown OpenCASCADE point with tag %d", centerTag);
923 return false;
924 }
925 if(!_tagVertex.IsBound(endTag)) {
926 Msg::Error("Unknown OpenCASCADE point with tag %d", endTag);
927 return false;
928 }
929
930 TopoDS_Edge result;
931 try {
932 TopoDS_Vertex start = TopoDS::Vertex(_tagVertex.Find(startTag));
933 TopoDS_Vertex center = TopoDS::Vertex(_tagVertex.Find(centerTag));
934 TopoDS_Vertex end = TopoDS::Vertex(_tagVertex.Find(endTag));
935 gp_Pnt aP1 = BRep_Tool::Pnt(start);
936 gp_Pnt aP2 = BRep_Tool::Pnt(center);
937 gp_Pnt aP3 = BRep_Tool::Pnt(end);
938 Standard_Real Radius = aP1.Distance(aP2);
939 gce_MakeCirc MC(aP2, gce_MakePln(aP1, aP2, aP3).Value(), Radius);
940 if(!MC.IsDone()) {
941 Msg::Error("Could not build circle");
942 return false;
943 }
944 const gp_Circ &Circ = MC.Value();
945 Standard_Real Alpha1 = ElCLib::Parameter(Circ, aP1);
946 Standard_Real Alpha2 = ElCLib::Parameter(Circ, aP3);
947 Handle(Geom_Circle) C = new Geom_Circle(Circ);
948 Handle(Geom_TrimmedCurve) arc =
949 new Geom_TrimmedCurve(C, Alpha1, Alpha2, false);
950 BRepBuilderAPI_MakeEdge e(arc, start, end);
951 e.Build();
952 if(!e.IsDone()) {
953 Msg::Error("Could not create circle arc");
954 return false;
955 }
956 result = e.Edge();
957 } catch(Standard_Failure &err) {
958 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
959 return false;
960 }
961 if(tag < 0) tag = getMaxTag(1) + 1;
962 _bind(result, tag, true);
963 return true;
964 }
965
addEllipseArc(int & tag,int startTag,int centerTag,int majorTag,int endTag)966 bool OCC_Internals::addEllipseArc(int &tag, int startTag, int centerTag,
967 int majorTag, int endTag)
968 {
969 if(tag >= 0 && _tagEdge.IsBound(tag)) {
970 Msg::Error("OpenCASCADE curve with tag %d already exists", tag);
971 return false;
972 }
973 if(!_tagVertex.IsBound(startTag)) {
974 Msg::Error("Unknown OpenCASCADE point with tag %d", startTag);
975 return false;
976 }
977 if(!_tagVertex.IsBound(centerTag)) {
978 Msg::Error("Unknown OpenCASCADE point with tag %d", centerTag);
979 return false;
980 }
981 if(!_tagVertex.IsBound(majorTag)) {
982 Msg::Error("Unknown OpenCASCADE point with tag %d", majorTag);
983 return false;
984 }
985 if(!_tagVertex.IsBound(endTag)) {
986 Msg::Error("Unknown OpenCASCADE point with tag %d", endTag);
987 return false;
988 }
989
990 TopoDS_Edge result;
991 try {
992 TopoDS_Vertex start = TopoDS::Vertex(_tagVertex.Find(startTag));
993 TopoDS_Vertex center = TopoDS::Vertex(_tagVertex.Find(centerTag));
994 TopoDS_Vertex major = TopoDS::Vertex(_tagVertex.Find(majorTag));
995 TopoDS_Vertex end = TopoDS::Vertex(_tagVertex.Find(endTag));
996 gp_Pnt startPnt = BRep_Tool::Pnt(start);
997 gp_Pnt centerPnt = BRep_Tool::Pnt(center);
998 gp_Pnt majorPnt = BRep_Tool::Pnt(major);
999 gp_Pnt endPnt = BRep_Tool::Pnt(end);
1000 gp_XYZ x1 = startPnt.XYZ() - centerPnt.XYZ();
1001 gp_XYZ x2 = endPnt.XYZ() - centerPnt.XYZ();
1002 gp_Dir u = majorPnt.XYZ() - centerPnt.XYZ();
1003 gp_Dir v;
1004 if(!u.IsParallel(x1, 1e-6))
1005 v = x1 - x1.Dot(u.XYZ()) * u.XYZ();
1006 else if(!u.IsParallel(x2, 1e-6))
1007 v = x2 - x2.Dot(u.XYZ()) * u.XYZ();
1008 else {
1009 Msg::Error("The points do not define an ellipse");
1010 return false;
1011 }
1012 Standard_Real x1u = Square(x1.Dot(u.XYZ()));
1013 Standard_Real x1v = Square(x1.Dot(v.XYZ()));
1014 Standard_Real x2u = Square(x2.Dot(u.XYZ()));
1015 Standard_Real x2v = Square(x2.Dot(v.XYZ()));
1016 if(IsEqual(x1u, x2u) || IsEqual(x1v, x2v)) {
1017 Msg::Error("The points do not define an ellipse");
1018 return false;
1019 }
1020 Standard_Real a2 = (x1v * x2u - x1u * x2v) / (x1v - x2v);
1021 Standard_Real b2 = (x1u * x2v - x1v * x2u) / (x1u - x2u);
1022 if(a2 <= 0.0 || b2 <= 0.0) {
1023 Msg::Error("The points do not define an ellipse");
1024 return false;
1025 }
1026 Standard_Real a; // Major radius
1027 Standard_Real b; // Minor radius
1028 gp_Ax2 Axes; // Ellipse local coordinate system
1029 if(a2 >= b2) {
1030 a = Sqrt(a2);
1031 b = Sqrt(b2);
1032 Axes = gp_Ax2(centerPnt, u ^ v, u);
1033 }
1034 else {
1035 Msg::Warning("Major radius smaller than minor radius");
1036 a = Sqrt(b2);
1037 b = Sqrt(a2);
1038 Axes = gp_Ax2(centerPnt, v ^ u, v);
1039 }
1040 gce_MakeElips ME(Axes, a, b);
1041 if(!ME.IsDone()) {
1042 Msg::Error("Could not build ellipse");
1043 return false;
1044 }
1045 const gp_Elips &Elips = ME.Value();
1046 Standard_Real Alpha1 = ElCLib::Parameter(Elips, startPnt);
1047 Standard_Real Alpha2 = ElCLib::Parameter(Elips, endPnt);
1048 Handle(Geom_Ellipse) E = new Geom_Ellipse(Elips);
1049 Handle(Geom_TrimmedCurve) arc;
1050 if((Alpha2 > Alpha1 && Alpha2 - Alpha1 < M_PI) || Alpha1 - Alpha2 > M_PI)
1051 arc = new Geom_TrimmedCurve(E, Alpha1, Alpha2, true);
1052 else
1053 arc = new Geom_TrimmedCurve(E, Alpha2, Alpha1, false);
1054 BRepBuilderAPI_MakeEdge e(arc, start, end);
1055 e.Build();
1056 if(!e.IsDone()) {
1057 Msg::Error("Could not create ellipse arc");
1058 return false;
1059 }
1060 result = e.Edge();
1061 } catch(Standard_Failure &err) {
1062 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1063 return false;
1064 }
1065 if(tag < 0) tag = getMaxTag(1) + 1;
1066 _bind(result, tag, true);
1067 return true;
1068 }
1069
addCircle(int & tag,double x,double y,double z,double r,double angle1,double angle2)1070 bool OCC_Internals::addCircle(int &tag, double x, double y, double z, double r,
1071 double angle1, double angle2)
1072 {
1073 if(tag >= 0 && _tagEdge.IsBound(tag)) {
1074 Msg::Error("OpenCASCADE curve with tag %d already exists", tag);
1075 return false;
1076 }
1077 if(r <= 0) {
1078 Msg::Error("Circle radius should be positive");
1079 return false;
1080 }
1081
1082 TopoDS_Edge result;
1083 try {
1084 gp_Dir N_dir(0., 0., 1.);
1085 gp_Dir x_dir(1., 0., 0.);
1086 gp_Pnt center(x, y, z);
1087 gp_Ax2 axis(center, N_dir, x_dir);
1088 gp_Circ circ(axis, r);
1089 if(angle1 == 0. && angle2 == 2 * M_PI) {
1090 result = BRepBuilderAPI_MakeEdge(circ);
1091 }
1092 else {
1093 Handle(Geom_Circle) C = new Geom_Circle(circ);
1094 Handle(Geom_TrimmedCurve) arc =
1095 new Geom_TrimmedCurve(C, angle1, angle2, false);
1096 BRepBuilderAPI_MakeEdge e(arc);
1097 if(!e.IsDone()) {
1098 Msg::Error("Could not create circle arc");
1099 return false;
1100 }
1101 result = e.Edge();
1102 }
1103 } catch(Standard_Failure &err) {
1104 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1105 return false;
1106 }
1107 if(tag < 0) tag = getMaxTag(1) + 1;
1108 _bind(result, tag, true);
1109 return true;
1110 }
1111
addEllipse(int & tag,double x,double y,double z,double rx,double ry,double angle1,double angle2)1112 bool OCC_Internals::addEllipse(int &tag, double x, double y, double z,
1113 double rx, double ry, double angle1,
1114 double angle2)
1115 {
1116 if(tag >= 0 && _tagEdge.IsBound(tag)) {
1117 Msg::Error("OpenCASCADE curve with tag %d already exists", tag);
1118 return false;
1119 }
1120 if(ry > rx) {
1121 Msg::Error("Major radius rx should be larger than minor radius ry");
1122 return false;
1123 }
1124 if(ry <= 0 || rx <= 0) {
1125 Msg::Error("Ellipse radii should be positive");
1126 return false;
1127 }
1128
1129 TopoDS_Edge result;
1130 try {
1131 gp_Dir N_dir(0., 0., 1.);
1132 gp_Dir x_dir(1., 0., 0.);
1133 gp_Pnt center(x, y, z);
1134 gp_Ax2 axis(center, N_dir, x_dir);
1135 gp_Elips elips(axis, rx, ry);
1136 if(angle1 == 0 && angle2 == 2 * M_PI) {
1137 result = BRepBuilderAPI_MakeEdge(elips);
1138 }
1139 else {
1140 Handle(Geom_Ellipse) E = new Geom_Ellipse(elips);
1141 Handle(Geom_TrimmedCurve) arc =
1142 new Geom_TrimmedCurve(E, angle1, angle2, true);
1143 BRepBuilderAPI_MakeEdge e(arc);
1144 if(!e.IsDone()) {
1145 Msg::Error("Could not create ellipse arc");
1146 return false;
1147 }
1148 result = e.Edge();
1149 }
1150 } catch(Standard_Failure &err) {
1151 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1152 return false;
1153 }
1154 if(tag < 0) tag = getMaxTag(1) + 1;
1155 _bind(result, tag, true);
1156 return true;
1157 }
1158
debugBSpline(const Handle (Geom_BSplineCurve)& curve)1159 void debugBSpline(const Handle(Geom_BSplineCurve) & curve)
1160 {
1161 int degree = curve->Degree();
1162 bool periodic = curve->IsPeriodic();
1163 bool rational = curve->IsRational();
1164
1165 int npoles = curve->NbPoles();
1166 TColgp_Array1OfPnt poles(1, npoles);
1167 curve->Poles(poles);
1168
1169 TColStd_Array1OfReal weights(1, npoles);
1170 curve->Weights(weights);
1171
1172 int nknots = curve->NbKnots();
1173 TColStd_Array1OfReal knots(1, nknots);
1174 curve->Knots(knots);
1175
1176 TColStd_Array1OfInteger mults(1, nknots);
1177 curve->Multiplicities(mults);
1178
1179 printf("BSpline: degree %d, periodic %d, rational %d\n", degree, periodic,
1180 rational);
1181 printf("Poles:\n");
1182 for(int i = 1; i <= npoles; i++)
1183 printf(" %d (%g, %g, %g) weight %g\n", i, poles(i).X(), poles(i).Y(),
1184 poles(i).Z(), weights(i));
1185 printf("Knots:\n");
1186 for(int i = 1; i <= nknots; i++)
1187 printf(" %d (%g) mult %d\n", i, knots(i), mults(i));
1188 }
1189
_addBSpline(int & tag,const std::vector<int> & pointTags,int mode,const int degree,const std::vector<double> & weights,const std::vector<double> & knots,const std::vector<int> & multiplicities)1190 bool OCC_Internals::_addBSpline(int &tag, const std::vector<int> &pointTags,
1191 int mode, const int degree,
1192 const std::vector<double> &weights,
1193 const std::vector<double> &knots,
1194 const std::vector<int> &multiplicities)
1195 {
1196 if(tag >= 0 && _tagEdge.IsBound(tag)) {
1197 Msg::Error("OpenCASCADE curve with tag %d already exists", tag);
1198 return false;
1199 }
1200 if(pointTags.size() < 2) {
1201 Msg::Error("Number of control points should be at least 2");
1202 return false;
1203 }
1204
1205 TopoDS_Edge result;
1206 try {
1207 TColgp_Array1OfPnt ctrlPoints(1, pointTags.size());
1208 TopoDS_Vertex start, end;
1209 for(std::size_t i = 0; i < pointTags.size(); i++) {
1210 if(!_tagVertex.IsBound(pointTags[i])) {
1211 Msg::Error("Unknown OpenCASCADE point with tag %d", pointTags[i]);
1212 return false;
1213 }
1214 TopoDS_Vertex vertex = TopoDS::Vertex(_tagVertex.Find(pointTags[i]));
1215 ctrlPoints.SetValue(i + 1, BRep_Tool::Pnt(vertex));
1216 if(i == 0) start = vertex;
1217 if(i == pointTags.size() - 1) end = vertex;
1218 }
1219 bool periodic = (pointTags.front() == pointTags.back());
1220 if(mode == 0) {
1221 // BSpline through points (called "Spline" in Gmsh; will be C2, whereas it
1222 // is only C1 in the GEO kernel)
1223 int np = periodic ? ctrlPoints.Length() - 1 : ctrlPoints.Length();
1224 Handle(TColgp_HArray1OfPnt) p = new TColgp_HArray1OfPnt(1, np);
1225 for(int i = 1; i <= np; i++) p->SetValue(i, ctrlPoints(i));
1226 GeomAPI_Interpolate intp(p, periodic, CTX::instance()->geom.tolerance);
1227 intp.Perform();
1228 if(!intp.IsDone()) {
1229 Msg::Error("Could not interpolate spline");
1230 return false;
1231 }
1232 Handle(Geom_BSplineCurve) curve = intp.Curve();
1233 BRepBuilderAPI_MakeEdge e(curve, start, end);
1234 if(!e.IsDone()) {
1235 Msg::Error("Could not create spline");
1236 return false;
1237 }
1238 result = e.Edge();
1239 }
1240 else if(mode == 1) {
1241 // Bezier curve
1242 Handle(Geom_BezierCurve) curve = new Geom_BezierCurve(ctrlPoints);
1243 BRepBuilderAPI_MakeEdge e(curve, start, end);
1244 if(!e.IsDone()) {
1245 Msg::Error("Could not create Bezier curve");
1246 return false;
1247 }
1248 result = e.Edge();
1249 }
1250 else if(mode == 2) {
1251 // General BSpline curve, polynomial or rational, with explicit degree,
1252 // weights, knots and multiplicities
1253 if(degree < 0) {
1254 Msg::Error("BSpline degree (%d) should be >= 0", degree);
1255 return false;
1256 }
1257 if(weights.size() != pointTags.size()) {
1258 Msg::Error("Number of BSpline weights (%d) and control points (%d) "
1259 "should be equal",
1260 weights.size(), pointTags.size());
1261 return false;
1262 }
1263 if(knots.size() != multiplicities.size()) {
1264 Msg::Error(
1265 "Number of BSpline knots (%d) and multiplicities (%d) should be "
1266 "equal",
1267 knots.size(), multiplicities.size());
1268 return false;
1269 }
1270 if(knots.size() < 2) {
1271 Msg::Error("Number of BSpline knots (%d) should be >= 2", knots.size());
1272 return false;
1273 }
1274 for(std::size_t i = 0; i < knots.size() - 1; i++) {
1275 if(knots[i] >= knots[i + 1]) {
1276 Msg::Error("BSpline knots should be increasing: knot %d (%g) > "
1277 "knot %d (%g)",
1278 i, knots[i], i + 1, knots[i + 1]);
1279 return false;
1280 }
1281 }
1282 for(std::size_t i = 0; i < multiplicities.size(); i++) {
1283 if(multiplicities[i] < 1) {
1284 Msg::Error("BSpline multiplicities should be >= 1");
1285 return false;
1286 }
1287 if(i != 0 && i != multiplicities.size() - 1 &&
1288 multiplicities[i] > degree) {
1289 Msg::Error(
1290 "BSpline interior knot multiplicities should be <= degree");
1291 return false;
1292 }
1293 if((i == 0 || i == multiplicities.size() - 1) &&
1294 multiplicities[i] > degree + 1) {
1295 Msg::Error("BSpline end knot multiplicities should be <= degree + 1");
1296 return false;
1297 }
1298 }
1299 if(periodic) {
1300 if(multiplicities.front() != multiplicities.back()) {
1301 Msg::Error(
1302 "Periodic BSpline end knot multiplicies (%d and %d) should "
1303 "be equal",
1304 multiplicities.front(), multiplicities.back());
1305 return false;
1306 }
1307 const auto sum = std::accumulate(
1308 begin(multiplicities), std::prev(std::end(multiplicities)),
1309 std::size_t(0),
1310 [](std::size_t const partial_sum, int const multiplicity) {
1311 return partial_sum + multiplicity;
1312 });
1313 if(pointTags.size() - 1 != sum) {
1314 Msg::Error("Number of control points - 1 for periodic BSpline should "
1315 "be equal to the sum of multiplicities for all knots "
1316 "except the first (or last)");
1317 return false;
1318 }
1319 }
1320 else {
1321 const auto sum = std::accumulate(
1322 begin(multiplicities), std::end(multiplicities), std::size_t(0),
1323 [](std::size_t const partial_sum, int const multiplicity) {
1324 return partial_sum + multiplicity;
1325 });
1326 if(pointTags.size() != sum - degree - 1) {
1327 Msg::Error("Number of control points for non-periodic BSpline should "
1328 "be equal to the sum of multiplicities - degree - 1");
1329 return false;
1330 }
1331 }
1332 int np = (periodic ? ctrlPoints.Length() - 1 : ctrlPoints.Length());
1333 TColgp_Array1OfPnt p(1, np);
1334 TColStd_Array1OfReal w(1, np);
1335 for(int i = 1; i <= np; i++) {
1336 p.SetValue(i, ctrlPoints(i));
1337 w.SetValue(i, weights[i - 1]);
1338 }
1339 TColStd_Array1OfReal k(1, knots.size());
1340 for(std::size_t i = 0; i < knots.size(); i++) k.SetValue(i + 1, knots[i]);
1341 TColStd_Array1OfInteger m(1, multiplicities.size());
1342 for(std::size_t i = 0; i < multiplicities.size(); i++)
1343 m.SetValue(i + 1, multiplicities[i]);
1344 Handle(Geom_BSplineCurve) curve =
1345 new Geom_BSplineCurve(p, w, k, m, degree, periodic);
1346 if(curve->StartPoint().IsEqual(BRep_Tool::Pnt(start),
1347 CTX::instance()->geom.tolerance) &&
1348 curve->EndPoint().IsEqual(BRep_Tool::Pnt(end),
1349 CTX::instance()->geom.tolerance)) {
1350 BRepBuilderAPI_MakeEdge e(curve, start, end);
1351 if(!e.IsDone()) {
1352 Msg::Error("Could not create BSpline curve (with end points)");
1353 return false;
1354 }
1355 result = e.Edge();
1356 }
1357 else { // will create new topo vertices as necessary
1358 BRepBuilderAPI_MakeEdge e(curve);
1359 if(!e.IsDone()) {
1360 Msg::Error("Could not create BSpline curve (without end points)");
1361 return false;
1362 }
1363 result = e.Edge();
1364 // copy mesh size from start control point to new topo vertex; this way
1365 // we'll behave more like the built-in kernel (although the built-in
1366 // kernel keeps track of all the control points)
1367 double lc = _attributes->getMeshSize(0, start);
1368 if(lc > 0 && lc < MAX_LC) {
1369 TopExp_Explorer exp0;
1370 for(exp0.Init(result, TopAbs_VERTEX); exp0.More(); exp0.Next()) {
1371 TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
1372 _attributes->insert(new OCCAttributes(0, vertex, lc));
1373 }
1374 }
1375 }
1376 }
1377 } catch(Standard_Failure &err) {
1378 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1379 return false;
1380 }
1381 if(tag < 0) tag = getMaxTag(1) + 1;
1382 _bind(result, tag, true);
1383 return true;
1384 }
1385
addSpline(int & tag,const std::vector<int> & pointTags)1386 bool OCC_Internals::addSpline(int &tag, const std::vector<int> &pointTags)
1387 {
1388 return _addBSpline(tag, pointTags, 0);
1389 }
1390
addBezier(int & tag,const std::vector<int> & pointTags)1391 bool OCC_Internals::addBezier(int &tag, const std::vector<int> &pointTags)
1392 {
1393 return _addBSpline(tag, pointTags, 1);
1394 }
1395
addBSpline(int & tag,const std::vector<int> & pointTags,const int degree,const std::vector<double> & weights,const std::vector<double> & knots,const std::vector<int> & multiplicities)1396 bool OCC_Internals::addBSpline(int &tag, const std::vector<int> &pointTags,
1397 const int degree,
1398 const std::vector<double> &weights,
1399 const std::vector<double> &knots,
1400 const std::vector<int> &multiplicities)
1401 {
1402 int np = pointTags.size();
1403 if(np < 2) {
1404 Msg::Error("BSpline curve requires at least 2 control points");
1405 return false;
1406 }
1407 int d = degree;
1408 std::vector<double> w(weights), k(knots);
1409 std::vector<int> m(multiplicities);
1410 // degree 3 if not specified...
1411 if(d <= 0) d = 3;
1412 // ... or number of control points - 1 if not enough points
1413 if(d > np - 1) d = np - 1;
1414 // automatic default weights if not provided:
1415 if(w.empty()) w.resize(np, 1);
1416 // automatic default knots and multiplicities if not provided:
1417 if(k.empty()) {
1418 bool periodic = (pointTags.front() == pointTags.back());
1419 if(!periodic) {
1420 int sumOfAllMult = np + d + 1;
1421 int numKnots = sumOfAllMult - 2 * d;
1422 if(numKnots < 2) {
1423 Msg::Error("Not enough control points for building BSpline of "
1424 "degree %d",
1425 d);
1426 return false;
1427 }
1428 k.resize(numKnots);
1429 for(std::size_t i = 0; i < k.size(); i++) k[i] = i;
1430 m.resize(numKnots, 1);
1431 m.front() = d + 1;
1432 m.back() = d + 1;
1433 }
1434 else {
1435 k.resize(np - d + 2);
1436 for(std::size_t i = 0; i < k.size(); i++) k[i] = i;
1437 m.resize(k.size(), 1);
1438 m.front() = d - 1;
1439 m.back() = d - 1;
1440 }
1441 }
1442 return _addBSpline(tag, pointTags, 2, d, w, k, m);
1443 }
1444
addWire(int & tag,const std::vector<int> & curveTags,bool checkClosed)1445 bool OCC_Internals::addWire(int &tag, const std::vector<int> &curveTags,
1446 bool checkClosed)
1447 {
1448 if(tag >= 0 && _tagWire.IsBound(tag)) {
1449 Msg::Error("OpenCASCADE wire or curve loop with tag %d already exists", tag);
1450 return false;
1451 }
1452
1453 // Note: contrary to shells wires are always "sewed", i.e., a valid wire is
1454 // constructed if points are geometrically at the same location (even if they
1455 // are not topologically identical); there is thus no need to add a "sewing"
1456 // option.
1457 try {
1458 BRepBuilderAPI_MakeWire w;
1459 TopoDS_Wire wire;
1460 for(std::size_t i = 0; i < curveTags.size(); i++) {
1461 // all curve tags are > 0 for OCC : but to improve compatibility between
1462 // GEO and OCC factories, we allow negative tags - and simply ignore the
1463 // sign here
1464 int t = std::abs(curveTags[i]);
1465 if(!_tagEdge.IsBound(t)) {
1466 Msg::Error("Unknown OpenCASCADE curve with tag %d", t);
1467 return false;
1468 }
1469 TopoDS_Edge edge = TopoDS::Edge(_tagEdge.Find(t));
1470 w.Add(edge);
1471 }
1472 wire = w.Wire();
1473 if(checkClosed && !wire.Closed()) {
1474 Msg::Error("Curve loop is not closed");
1475 return false;
1476 }
1477 if(tag < 0) tag = getMaxTag(-1) + 1;
1478 _bind(wire, tag, true);
1479 } catch(Standard_Failure &err) {
1480 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1481 return false;
1482 }
1483 return true;
1484 }
1485
addCurveLoop(int & tag,const std::vector<int> & curveTags)1486 bool OCC_Internals::addCurveLoop(int &tag, const std::vector<int> &curveTags)
1487 {
1488 return addWire(tag, curveTags, true);
1489 }
1490
makeRectangle(TopoDS_Face & result,double x,double y,double z,double dx,double dy,double roundedRadius)1491 static bool makeRectangle(TopoDS_Face &result, double x, double y, double z,
1492 double dx, double dy, double roundedRadius)
1493 {
1494 if(!dx || !dy) {
1495 Msg::Error("Rectangle with zero width or height");
1496 return false;
1497 }
1498 try {
1499 TopoDS_Wire wire;
1500 if(roundedRadius <= 0.) {
1501 double x1 = x, y1 = y, z1 = z, x2 = x1 + dx, y2 = y1 + dy;
1502 TopoDS_Vertex v1 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1, y1, z1));
1503 TopoDS_Vertex v2 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2, y1, z1));
1504 TopoDS_Vertex v3 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2, y2, z1));
1505 TopoDS_Vertex v4 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1, y2, z1));
1506 TopoDS_Edge e1 = BRepBuilderAPI_MakeEdge(v1, v2);
1507 TopoDS_Edge e2 = BRepBuilderAPI_MakeEdge(v2, v3);
1508 TopoDS_Edge e3 = BRepBuilderAPI_MakeEdge(v3, v4);
1509 TopoDS_Edge e4 = BRepBuilderAPI_MakeEdge(v4, v1);
1510 wire = BRepBuilderAPI_MakeWire(e1, e2, e3, e4);
1511 }
1512 else {
1513 double x1, y1, z1 = z, x2, y2;
1514 double r = roundedRadius;
1515 if(dx > 0.) {
1516 x1 = x;
1517 x2 = x1 + dx;
1518 }
1519 else {
1520 x2 = x;
1521 x1 = x2 + dx;
1522 }
1523 if(dy > 0.) {
1524 y1 = y;
1525 y2 = y1 + dy;
1526 }
1527 else {
1528 y2 = y;
1529 y1 = y2 + dy;
1530 }
1531 TopoDS_Vertex v1 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1 + r, y1, z1));
1532 TopoDS_Vertex v2 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2 - r, y1, z1));
1533 TopoDS_Vertex v3 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2, y1 + r, z1));
1534 TopoDS_Vertex v4 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2, y2 - r, z1));
1535 TopoDS_Vertex v5 = BRepBuilderAPI_MakeVertex(gp_Pnt(x2 - r, y2, z1));
1536 TopoDS_Vertex v6 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1 + r, y2, z1));
1537 TopoDS_Vertex v7 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1, y2 - r, z1));
1538 TopoDS_Vertex v8 = BRepBuilderAPI_MakeVertex(gp_Pnt(x1, y1 + r, z1));
1539 TopoDS_Edge e1 = BRepBuilderAPI_MakeEdge(v1, v2);
1540 TopoDS_Edge e2 = BRepBuilderAPI_MakeEdge(v3, v4);
1541 TopoDS_Edge e3 = BRepBuilderAPI_MakeEdge(v5, v6);
1542 TopoDS_Edge e4 = BRepBuilderAPI_MakeEdge(v7, v8);
1543 gp_Pnt c1(x1 + r, y1 + r, z1);
1544 gp_Pnt c2(x2 - r, y1 + r, z1);
1545 gp_Pnt c3(x2 - r, y2 - r, z1);
1546 gp_Pnt c4(x1 + r, y2 - r, z1);
1547 gp_Pln plane = gce_MakePln(c1, c2, c3).Value();
1548 gp_Circ circ1 = gce_MakeCirc(c1, plane, r).Value();
1549 gp_Circ circ2 = gce_MakeCirc(c2, plane, r).Value();
1550 gp_Circ circ3 = gce_MakeCirc(c3, plane, r).Value();
1551 gp_Circ circ4 = gce_MakeCirc(c4, plane, r).Value();
1552 Handle(Geom_Circle) circle1 = new Geom_Circle(circ1);
1553 Handle(Geom_Circle) circle2 = new Geom_Circle(circ2);
1554 Handle(Geom_Circle) circle3 = new Geom_Circle(circ3);
1555 Handle(Geom_Circle) circle4 = new Geom_Circle(circ4);
1556 Handle(Geom_TrimmedCurve) arc1 =
1557 new Geom_TrimmedCurve(circle1, -M_PI, -M_PI / 2., true);
1558 Handle(Geom_TrimmedCurve) arc2 =
1559 new Geom_TrimmedCurve(circle2, -M_PI / 2, 0, true);
1560 Handle(Geom_TrimmedCurve) arc3 =
1561 new Geom_TrimmedCurve(circle3, 0, M_PI / 2, true);
1562 Handle(Geom_TrimmedCurve) arc4 =
1563 new Geom_TrimmedCurve(circle4, M_PI / 2, M_PI, true);
1564 TopoDS_Edge ce1 = BRepBuilderAPI_MakeEdge(arc1, v8, v1);
1565 TopoDS_Edge ce2 = BRepBuilderAPI_MakeEdge(arc2, v2, v3);
1566 TopoDS_Edge ce3 = BRepBuilderAPI_MakeEdge(arc3, v4, v5);
1567 TopoDS_Edge ce4 = BRepBuilderAPI_MakeEdge(arc4, v6, v7);
1568 BRepBuilderAPI_MakeWire w;
1569 w.Add(e1);
1570 w.Add(ce2);
1571 w.Add(e2);
1572 w.Add(ce3);
1573 w.Add(e3);
1574 w.Add(ce4);
1575 w.Add(e4);
1576 w.Add(ce1);
1577 wire = w.Wire();
1578 }
1579 result = BRepBuilderAPI_MakeFace(wire);
1580 } catch(Standard_Failure &err) {
1581 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1582 return false;
1583 }
1584 return true;
1585 }
1586
addRectangle(int & tag,double x,double y,double z,double dx,double dy,double roundedRadius)1587 bool OCC_Internals::addRectangle(int &tag, double x, double y, double z,
1588 double dx, double dy, double roundedRadius)
1589 {
1590 if(tag >= 0 && _tagFace.IsBound(tag)) {
1591 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
1592 return false;
1593 }
1594 TopoDS_Face result;
1595 if(!makeRectangle(result, x, y, z, dx, dy, roundedRadius)) return false;
1596 if(tag < 0) tag = getMaxTag(2) + 1;
1597 _bind(result, tag, true);
1598 return true;
1599 }
1600
makeDisk(TopoDS_Face & result,double xc,double yc,double zc,double rx,double ry)1601 static bool makeDisk(TopoDS_Face &result, double xc, double yc, double zc,
1602 double rx, double ry)
1603 {
1604 if(ry > rx) {
1605 Msg::Error("Major radius rx should be larger than minor radius ry");
1606 return false;
1607 }
1608 if(ry <= 0 || rx <= 0) {
1609 Msg::Error("Disk radius should be positive");
1610 return false;
1611 }
1612 try {
1613 gp_Dir N_dir(0., 0., 1.);
1614 gp_Dir x_dir(1., 0., 0.);
1615 gp_Pnt center(xc, yc, zc);
1616 gp_Ax2 axis(center, N_dir, x_dir);
1617 gp_Elips ellipse = gp_Elips(axis, rx, ry);
1618 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge(ellipse);
1619 TopoDS_Wire wire = BRepBuilderAPI_MakeWire(edge);
1620 result = BRepBuilderAPI_MakeFace(wire);
1621 } catch(Standard_Failure &err) {
1622 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1623 return false;
1624 }
1625 return true;
1626 }
1627
addDisk(int & tag,double xc,double yc,double zc,double rx,double ry)1628 bool OCC_Internals::addDisk(int &tag, double xc, double yc, double zc,
1629 double rx, double ry)
1630 {
1631 if(tag >= 0 && _tagFace.IsBound(tag)) {
1632 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
1633 return false;
1634 }
1635 TopoDS_Face result;
1636 if(!makeDisk(result, xc, yc, zc, rx, ry)) return false;
1637 if(tag < 0) tag = getMaxTag(2) + 1;
1638 _bind(result, tag, true);
1639 return true;
1640 }
1641
addPlaneSurface(int & tag,const std::vector<int> & wireTags)1642 bool OCC_Internals::addPlaneSurface(int &tag, const std::vector<int> &wireTags)
1643 {
1644 if(tag >= 0 && _tagFace.IsBound(tag)) {
1645 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
1646 return false;
1647 }
1648
1649 std::vector<TopoDS_Wire> wires;
1650 for(std::size_t i = 0; i < wireTags.size(); i++) {
1651 // all wire tags are > 0 for OCC : to improve compatibility between GEO and
1652 // OCC factories, allow negative tags - and simply ignore the sign here
1653 int wireTag = std::abs(wireTags[i]);
1654 if(!_tagWire.IsBound(wireTag)) {
1655 Msg::Error("Unknown OpenCASCADE curve loop with tag %d", wireTag);
1656 return false;
1657 }
1658 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTag));
1659 wires.push_back(wire);
1660 }
1661
1662 TopoDS_Face result;
1663 if(wires.size() == 0) {
1664 Msg::Error("Plane surface requires at least one curve loop");
1665 return false;
1666 }
1667
1668 try {
1669 BRepBuilderAPI_MakeFace f(wires[0]);
1670 for(std::size_t i = 1; i < wires.size(); i++) {
1671 // holes
1672 TopoDS_Wire w = wires[i];
1673 w.Orientation(TopAbs_REVERSED);
1674 f.Add(w);
1675 }
1676 f.Build();
1677 if(!f.IsDone()) {
1678 Msg::Error("Could not create surface");
1679 return false;
1680 }
1681 result = f.Face();
1682 if(CTX::instance()->geom.occAutoFix) {
1683 // make sure wires are oriented correctly
1684 ShapeFix_Face fix(result);
1685 fix.Perform();
1686 result = fix.Face();
1687 }
1688 } catch(Standard_Failure &err) {
1689 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1690 return false;
1691 }
1692
1693 if(tag < 0) tag = getMaxTag(2) + 1;
1694 _bind(result, tag, true);
1695 return true;
1696 }
1697
addSurfaceFilling(int & tag,int wireTag,const std::vector<int> & pointTags,const std::vector<int> & surfaceTags,const std::vector<int> & surfaceContinuity,const int degree,const int numPointsOnCurves,const int numIter,const bool anisotropic,const double tol2d,const double tol3d,const double tolAng,const double tolCurv,const int maxDegree,const int maxSegments)1698 bool OCC_Internals::addSurfaceFilling(int &tag, int wireTag,
1699 const std::vector<int> &pointTags,
1700 const std::vector<int> &surfaceTags,
1701 const std::vector<int> &surfaceContinuity,
1702 const int degree,
1703 const int numPointsOnCurves,
1704 const int numIter,
1705 const bool anisotropic,
1706 const double tol2d,
1707 const double tol3d,
1708 const double tolAng,
1709 const double tolCurv,
1710 const int maxDegree,
1711 const int maxSegments)
1712 {
1713 if(tag >= 0 && _tagFace.IsBound(tag)) {
1714 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
1715 return false;
1716 }
1717
1718 TopoDS_Face result;
1719 try {
1720 BRepOffsetAPI_MakeFilling f(degree, numPointsOnCurves, numIter, anisotropic,
1721 tol2d, tol3d, tolAng, tolCurv, maxDegree,
1722 maxSegments);
1723 // bounding edge constraints
1724 if(!_tagWire.IsBound(wireTag)) {
1725 Msg::Error("Unknown OpenCASCADE curve loop with tag %d", wireTag);
1726 return false;
1727 }
1728 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTag));
1729 TopExp_Explorer exp0;
1730 std::size_t i = 0;
1731 for(exp0.Init(wire, TopAbs_EDGE); exp0.More(); exp0.Next()) {
1732 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
1733 if(i < surfaceTags.size()) {
1734 // associated face constraint (does not seem to work...)
1735 if(!_tagFace.IsBound(surfaceTags[i])) {
1736 Msg::Error("Unknown OpenCASCADE surface with tag %d", surfaceTags[i]);
1737 return false;
1738 }
1739 TopoDS_Face face = TopoDS::Face(_tagFace.Find(surfaceTags[i]));
1740 if(i < surfaceContinuity.size() && surfaceContinuity[i] == 2)
1741 f.Add(edge, face, GeomAbs_G2);
1742 else
1743 f.Add(edge, face, GeomAbs_G1);
1744 }
1745 else {
1746 f.Add(edge, GeomAbs_C0);
1747 }
1748 i++;
1749 }
1750 // point constraints
1751 for(std::size_t i = 0; i < pointTags.size(); i++) {
1752 if(!_tagVertex.IsBound(pointTags[i])) {
1753 Msg::Error("Unknown OpenCASCADE point with tag %d", pointTags[i]);
1754 return false;
1755 }
1756 TopoDS_Vertex vertex = TopoDS::Vertex(_tagVertex.Find(pointTags[i]));
1757 f.Add(BRep_Tool::Pnt(vertex));
1758 }
1759 f.Build();
1760 if(!f.IsDone()) {
1761 Msg::Error("Could not build surface filling");
1762 return false;
1763 }
1764 // face filling duplicates the edges, so we need to go back to the
1765 // underlying surface, and remake a new face explicitly with the wire;
1766 // applying ShapeFix is mandatory (not sure why...)
1767 TopoDS_Face tmp = TopoDS::Face(f.Shape());
1768 Handle(Geom_Surface) s = BRep_Tool::Surface(tmp);
1769 result = BRepBuilderAPI_MakeFace(s, wire);
1770 ShapeFix_Face fix(result);
1771 fix.SetPrecision(CTX::instance()->geom.tolerance);
1772 fix.Perform();
1773 fix.FixOrientation(); // and I don't understand why this is necessary
1774 result = fix.Face();
1775 } catch(Standard_Failure &err) {
1776 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1777 return false;
1778 }
1779
1780 if(tag < 0) tag = getMaxTag(2) + 1;
1781 _bind(result, tag, true);
1782 return true;
1783 }
1784
addBSplineFilling(int & tag,int wireTag,const std::string & type)1785 bool OCC_Internals::addBSplineFilling(int &tag, int wireTag,
1786 const std::string &type)
1787 {
1788 if(tag >= 0 && _tagFace.IsBound(tag)) {
1789 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
1790 return false;
1791 }
1792
1793 TopoDS_Face result;
1794 try {
1795 GeomFill_BSplineCurves f;
1796 if(!_tagWire.IsBound(wireTag)) {
1797 Msg::Error("Unknown OpenCASCADE curve loop with tag %d", wireTag);
1798 return false;
1799 }
1800 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTag));
1801 TopExp_Explorer exp0;
1802 std::vector<Handle(Geom_BSplineCurve)> bsplines;
1803 for(exp0.Init(wire, TopAbs_EDGE); exp0.More(); exp0.Next()) {
1804 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
1805 double s0, s1;
1806 Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, s0, s1);
1807 if(curve->DynamicType() == STANDARD_TYPE(Geom_BSplineCurve)) {
1808 bsplines.push_back(Handle(Geom_BSplineCurve)::DownCast(curve));
1809 }
1810 else {
1811 Msg::Error("Bounding curve for BSpline filling should be a BSpline");
1812 }
1813 }
1814
1815 GeomFill_FillingStyle t;
1816 if(type == "Stretch")
1817 t = GeomFill_StretchStyle; // flattest patch
1818 else if(type == "Coons")
1819 t = GeomFill_CoonsStyle; // rounded with less depth than Curved
1820 else
1821 t = GeomFill_CurvedStyle; // most rounded patch
1822
1823 if(bsplines.size() == 4) {
1824 f.Init(bsplines[0], bsplines[1], bsplines[2], bsplines[3], t);
1825 }
1826 else if(bsplines.size() == 3) {
1827 f.Init(bsplines[0], bsplines[1], bsplines[2], t);
1828 }
1829 else if(bsplines.size() == 2) {
1830 f.Init(bsplines[0], bsplines[1], t);
1831 }
1832 else {
1833 Msg::Error(
1834 "BSpline filling requires between 2 and 4 boundary BSpline curves");
1835 return false;
1836 }
1837 const Handle(Geom_BSplineSurface) &surf = f.Surface();
1838 result = BRepBuilderAPI_MakeFace(surf, wire);
1839 ShapeFix_Face fix(result); // not sure why, but this is necessary
1840 fix.SetPrecision(CTX::instance()->geom.tolerance);
1841 fix.Perform();
1842 fix.FixOrientation();
1843 result = fix.Face();
1844 } catch(Standard_Failure &err) {
1845 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1846 return false;
1847 }
1848
1849 if(tag < 0) tag = getMaxTag(2) + 1;
1850 _bind(result, tag, true);
1851 return true;
1852 }
1853
addBezierFilling(int & tag,int wireTag,const std::string & type)1854 bool OCC_Internals::addBezierFilling(int &tag, int wireTag,
1855 const std::string &type)
1856 {
1857 if(tag >= 0 && _tagFace.IsBound(tag)) {
1858 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
1859 return false;
1860 }
1861
1862 TopoDS_Face result;
1863 try {
1864 GeomFill_BezierCurves f;
1865 if(!_tagWire.IsBound(wireTag)) {
1866 Msg::Error("Unknown OpenCASCADE curve loop with tag %d", wireTag);
1867 return false;
1868 }
1869 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTag));
1870 TopExp_Explorer exp0;
1871 std::vector<Handle(Geom_BezierCurve)> beziers;
1872 for(exp0.Init(wire, TopAbs_EDGE); exp0.More(); exp0.Next()) {
1873 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
1874 double s0, s1;
1875 Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, s0, s1);
1876 if(curve->DynamicType() == STANDARD_TYPE(Geom_BezierCurve)) {
1877 beziers.push_back(Handle(Geom_BezierCurve)::DownCast(curve));
1878 }
1879 else {
1880 Msg::Error(
1881 "Bounding curve for Bezier filling should be a Bezier curve");
1882 }
1883 }
1884
1885 GeomFill_FillingStyle t;
1886 if(type == "Stretch")
1887 t = GeomFill_StretchStyle; // flattest patch
1888 else if(type == "Coons")
1889 t = GeomFill_CoonsStyle; // rounded with less depth than Curved
1890 else
1891 t = GeomFill_CurvedStyle; // most rounded patch
1892
1893 if(beziers.size() == 4) {
1894 f.Init(beziers[0], beziers[1], beziers[2], beziers[3], t);
1895 }
1896 else if(beziers.size() == 3) {
1897 f.Init(beziers[0], beziers[1], beziers[2], t);
1898 }
1899 else if(beziers.size() == 2) {
1900 f.Init(beziers[0], beziers[1], t);
1901 }
1902 else {
1903 Msg::Error(
1904 "Bezier filling requires between 2 and 4 boundary Bezier curves");
1905 return false;
1906 }
1907 const Handle(Geom_BezierSurface) &surf = f.Surface();
1908 result = BRepBuilderAPI_MakeFace(surf, wire);
1909 ShapeFix_Face fix(result); // not sure why, but this is necessary
1910 fix.SetPrecision(CTX::instance()->geom.tolerance);
1911 fix.Perform();
1912 fix.FixOrientation();
1913 result = fix.Face();
1914 } catch(Standard_Failure &err) {
1915 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1916 return false;
1917 }
1918
1919 if(tag < 0) tag = getMaxTag(2) + 1;
1920 _bind(result, tag, true);
1921 return true;
1922 }
1923
makeEdgeOnSurface(const TopoDS_Edge & edge,const Handle (Geom_Surface)& surf,bool curve3D,TopoDS_Edge & edgeOnSurf)1924 static bool makeEdgeOnSurface(const TopoDS_Edge &edge,
1925 const Handle(Geom_Surface) &surf,
1926 bool curve3D, TopoDS_Edge &edgeOnSurf)
1927 {
1928 try {
1929 Standard_Real first, last;
1930 Handle(Geom_Curve) c = BRep_Tool::Curve(edge, first, last);
1931 if(curve3D) {
1932 // use the 3D curves in the wire and project them onto the patch
1933 Handle(Geom_Curve) cProj = GeomProjLib::Project
1934 (new Geom_TrimmedCurve(c, first, last, Standard_True, Standard_False),
1935 surf);
1936 edgeOnSurf = BRepBuilderAPI_MakeEdge
1937 (cProj, cProj->FirstParameter(), cProj->LastParameter());
1938 }
1939 else {
1940 // assume the 3D curve is actually a 2D curve in the parametric plane of the
1941 // surface: to retrieve the 2D curve, project the 3D curve on the z=0 plane
1942 Handle(Geom_Plane) p = new Geom_Plane(0, 0, 1, 0);
1943 TopLoc_Location loc;
1944 Handle(Geom2d_Curve) c2d =
1945 BRep_Tool::CurveOnSurface(edge, p, loc, first, last);
1946 // BRep_Tool::CurveOnPlane(edge, p, loc, first, last); // OCCT >= 7.2
1947 edgeOnSurf = BRepBuilderAPI_MakeEdge(c2d, surf, first, last);
1948 }
1949 }
1950 catch(Standard_Failure &err) {
1951 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
1952 return false;
1953 }
1954 return true;
1955 }
1956
makeTrimmedSurface(const Handle (Geom_Surface)& surf,const std::vector<TopoDS_Wire> & wires,bool wire3D,TopoDS_Face & result)1957 static bool makeTrimmedSurface(const Handle(Geom_Surface) &surf,
1958 const std::vector<TopoDS_Wire> &wires,
1959 bool wire3D, TopoDS_Face &result)
1960 {
1961 if(wires.empty()) { // natural bounds
1962 result = BRepBuilderAPI_MakeFace(surf, CTX::instance()->geom.tolerance);
1963 #if 0
1964 // Activate this to use input points as corners if they are on the corners
1965 // of the patch. (Since the natural "Replace(old_vertex, new_vertex)" on the
1966 // face does not work, we do it on each edge. Sigh...) Since when buiding
1967 // multi-patch models a fragment or sewing will eventually be necessary to
1968 // glue the patches, it's not that useful let's leave this commented out.
1969 ShapeBuild_ReShape rebuild;
1970 TopExp_Explorer exp0;
1971 for(exp0.Init(result, TopAbs_EDGE); exp0.More(); exp0.Next()) {
1972 TopoDS_Edge e = TopoDS::Edge(exp0.Current());
1973 TopoDS_Vertex v1 = TopExp::FirstVertex(e);
1974 TopoDS_Vertex v2 = TopExp::LastVertex(e);
1975 double s0, s1;
1976 Handle(Geom_Curve) curve = BRep_Tool::Curve(e, s0, s1);
1977 if(curve->DynamicType() == STANDARD_TYPE(Geom_BSplineCurve)){
1978 Handle(Geom_BSplineCurve) bs = Handle(Geom_BSplineCurve)::DownCast(curve);
1979 for(std::size_t i = 0; i < corners.size(); i++) {
1980 if(bs->StartPoint().IsEqual(BRep_Tool::Pnt(corners[i]),
1981 CTX::instance()->geom.tolerance)) {
1982 v1 = corners[i];
1983 }
1984 if(bs->EndPoint().IsEqual(BRep_Tool::Pnt(corners[i]),
1985 CTX::instance()->geom.tolerance)) {
1986 v2 = corners[i];
1987 }
1988 }
1989 }
1990 BRepBuilderAPI_MakeEdge newe(curve, v1, v2);
1991 rebuild.Replace(e, newe);
1992 }
1993 result = TopoDS::Face(rebuild.Apply(result));
1994 ShapeFix_Face fix(result); // not sure why, but this is necessary
1995 fix.SetPrecision(CTX::instance()->geom.tolerance);
1996 fix.Perform();
1997 fix.FixOrientation();
1998 result = fix.Face();
1999 #endif
2000 }
2001 else { // trimmed patch, using provided wires
2002 std::vector<TopoDS_Wire> wiresProj;
2003 for(std::size_t i = 0; i < wires.size(); i++) {
2004 BRepBuilderAPI_MakeWire w;
2005 TopExp_Explorer exp0;
2006 for(exp0.Init(wires[i], TopAbs_EDGE); exp0.More(); exp0.Next()) {
2007 TopoDS_Edge edge = TopoDS::Edge(exp0.Current()), edgeOnSurf;
2008 if(makeEdgeOnSurface(edge, surf, wire3D, edgeOnSurf))
2009 w.Add(edgeOnSurf);
2010 }
2011 TopoDS_Wire wire = w.Wire();
2012 wiresProj.push_back(wire);
2013 }
2014 BRepBuilderAPI_MakeFace f(surf, wiresProj[0]);
2015 for(std::size_t i = 1; i < wiresProj.size(); i++) f.Add(wiresProj[i]);
2016 result = f.Face();
2017 // recover 3D curves for pcurves
2018 ShapeFix_Face fix(result);
2019 fix.Perform();
2020 result = fix.Face();
2021 }
2022 return true;
2023 }
2024
addBSplineSurface(int & tag,const std::vector<int> & pointTags,const int numPointsU,const int degreeU,const int degreeV,const std::vector<double> & weights,const std::vector<double> & knotsU,const std::vector<double> & knotsV,const std::vector<int> & multiplicitiesU,const std::vector<int> & multiplicitiesV,const std::vector<int> & wireTags,bool wire3D)2025 bool OCC_Internals::addBSplineSurface(
2026 int &tag, const std::vector<int> &pointTags, const int numPointsU,
2027 const int degreeU, const int degreeV, const std::vector<double> &weights,
2028 const std::vector<double> &knotsU, const std::vector<double> &knotsV,
2029 const std::vector<int> &multiplicitiesU,
2030 const std::vector<int> &multiplicitiesV, const std::vector<int> &wireTags,
2031 bool wire3D)
2032 {
2033 if(tag >= 0 && _tagFace.IsBound(tag)) {
2034 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
2035 return false;
2036 }
2037
2038 // deal with default values
2039 if(numPointsU < 1) {
2040 Msg::Error("Wrong number of control points along U for BSpline surface");
2041 return false;
2042 }
2043 int numPointsV = pointTags.size() / numPointsU;
2044 if(numPointsU * numPointsV != (int)pointTags.size()) {
2045 Msg::Error("Wrong number of control points for BSpline surface");
2046 return false;
2047 }
2048 int dU = degreeU, dV = degreeV;
2049 std::vector<double> w(weights), kU(knotsU), kV(knotsV);
2050 std::vector<int> mU(multiplicitiesU), mV(multiplicitiesV);
2051 // degree 3 if not specified...
2052 if(dU <= 0) dU = 3;
2053 if(dV <= 0) dV = 3;
2054 // ... or number of control points - 1 if not enough points
2055 if(dU > numPointsU - 1) dU = numPointsU - 1;
2056 if(dV > numPointsV - 1) dV = numPointsV - 1;
2057 // automatic default weights if not provided:
2058 if(w.empty()) w.resize(pointTags.size(), 1);
2059 if(w.size() != pointTags.size()) {
2060 Msg::Error("Wrong number of weights for BSpline surface");
2061 return false;
2062 }
2063 bool periodicU = true;
2064 for(int i = 0; i < numPointsV; i++) {
2065 if(pointTags[i * numPointsU] != pointTags[(i + 1) * numPointsU - 1]) {
2066 periodicU = false;
2067 break;
2068 }
2069 }
2070 bool periodicV = true;
2071 for(int i = 0; i < numPointsU; i++) {
2072 if(pointTags[i * numPointsV] != pointTags[(i + 1) * numPointsV - 1]) {
2073 periodicV = false;
2074 break;
2075 }
2076 }
2077 // automatic default knots and multiplicities along U if not provided:
2078 if(kU.empty()) {
2079 if(!periodicU) {
2080 int sumOfAllMultU = numPointsU + dU + 1;
2081 int numKnotsU = sumOfAllMultU - 2 * dU;
2082 if(numKnotsU < 2) {
2083 Msg::Error("Not enough control points along U for building BSpline of "
2084 "degree %d x %d",
2085 dU, dV);
2086 return false;
2087 }
2088 kU.resize(numKnotsU);
2089 for(std::size_t i = 0; i < kU.size(); i++) kU[i] = i;
2090 mU.resize(numKnotsU, 1);
2091 mU.front() = dU + 1;
2092 mU.back() = dU + 1;
2093 }
2094 else {
2095 kU.resize(numPointsU - dU + 2);
2096 for(std::size_t i = 0; i < kU.size(); i++) kU[i] = i;
2097 mU.resize(kU.size(), 1);
2098 mU.front() = dU - 1;
2099 mU.back() = dU - 1;
2100 }
2101 }
2102 if(kU.size() != mU.size()) {
2103 Msg::Error("Number of BSpline knots and multiplicities should be equal");
2104 return false;
2105 }
2106 // automatic default knots and multiplicities along V if not provided:
2107 if(kV.empty()) {
2108 if(!periodicV) {
2109 int sumOfAllMultV = numPointsV + dV + 1;
2110 int numKnotsV = sumOfAllMultV - 2 * dV;
2111 if(numKnotsV < 2) {
2112 Msg::Error("Not enough control points along V for building BSpline of "
2113 "degree %d x %d",
2114 dU, dV);
2115 return false;
2116 }
2117 kV.resize(numKnotsV);
2118 for(std::size_t i = 0; i < kV.size(); i++) kV[i] = i;
2119 mV.resize(numKnotsV, 1);
2120 mV.front() = dV + 1;
2121 mV.back() = dV + 1;
2122 }
2123 else {
2124 kV.resize(numPointsV - dV + 2);
2125 for(std::size_t i = 0; i < kV.size(); i++) kV[i] = i;
2126 mV.resize(kV.size(), 1);
2127 mV.front() = dV - 1;
2128 mV.back() = dV - 1;
2129 }
2130 }
2131 if(kV.size() != mV.size()) {
2132 Msg::Error("Number of BSpline knots and multiplicities should be equal");
2133 return false;
2134 }
2135
2136 std::vector<TopoDS_Wire> wires;
2137 for(std::size_t i = 0; i < wireTags.size(); i++) {
2138 int wireTag = std::abs(wireTags[i]);
2139 if(!_tagWire.IsBound(wireTag)) {
2140 Msg::Error("Unknown OpenCASCADE curve loop with tag %d", wireTag);
2141 return false;
2142 }
2143 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTag));
2144 wires.push_back(wire);
2145 }
2146
2147 TopoDS_Face result;
2148 try {
2149 std::vector<TopoDS_Vertex> corners;
2150 int npU = (periodicU ? numPointsU - 1 : numPointsU);
2151 int npV = (periodicV ? numPointsV - 1 : numPointsV);
2152 TColgp_Array2OfPnt pp(1, npU, 1, npV);
2153 for(int i = 1; i <= npU; i++) {
2154 for(int j = 1; j <= npV; j++) {
2155 int k = (j - 1) * numPointsU + (i - 1);
2156 if(!_tagVertex.IsBound(pointTags[k])) {
2157 Msg::Error("Unknown OpenCASCADE point with tag %d", pointTags[k]);
2158 return false;
2159 }
2160 TopoDS_Vertex vertex = TopoDS::Vertex(_tagVertex.Find(pointTags[k]));
2161 if((i == 1 && j == 1) || (i == 1 && j == npV) || (i == npU && j == 1) ||
2162 (i == npU && j == npV))
2163 corners.push_back(vertex);
2164 pp.SetValue(i, j, BRep_Tool::Pnt(vertex));
2165 }
2166 }
2167 TColStd_Array2OfReal ww(1, npU, 1, npV);
2168 for(int i = 1; i <= npU; i++) {
2169 for(int j = 1; j <= npV; j++) {
2170 int k = (j - 1) * numPointsU + (i - 1);
2171 ww.SetValue(i, j, w[k]);
2172 }
2173 }
2174 TColStd_Array1OfReal kkU(1, kU.size());
2175 for(std::size_t i = 1; i <= kU.size(); i++) kkU.SetValue(i, kU[i - 1]);
2176 TColStd_Array1OfReal kkV(1, kV.size());
2177 for(std::size_t i = 1; i <= kV.size(); i++) kkV.SetValue(i, kV[i - 1]);
2178 TColStd_Array1OfInteger mmU(1, mU.size());
2179 for(std::size_t i = 1; i <= mU.size(); i++) mmU.SetValue(i, mU[i - 1]);
2180 TColStd_Array1OfInteger mmV(1, mV.size());
2181 for(std::size_t i = 1; i <= mV.size(); i++) mmV.SetValue(i, mV[i - 1]);
2182 Handle(Geom_BSplineSurface) surf = new Geom_BSplineSurface(
2183 pp, ww, kkU, kkV, mmU, mmV, dU, dV, periodicU, periodicV);
2184 makeTrimmedSurface(surf, wires, wire3D, result);
2185 } catch(Standard_Failure &err) {
2186 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2187 return false;
2188 }
2189
2190 if(tag < 0) tag = getMaxTag(2) + 1;
2191 _bind(result, tag, true);
2192 return true;
2193 }
2194
addBezierSurface(int & tag,const std::vector<int> & pointTags,const int numPointsU,const std::vector<int> & wireTags,bool wire3D)2195 bool OCC_Internals::addBezierSurface(int &tag,
2196 const std::vector<int> &pointTags,
2197 const int numPointsU,
2198 const std::vector<int> &wireTags,
2199 bool wire3D)
2200 {
2201 if(tag >= 0 && _tagFace.IsBound(tag)) {
2202 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
2203 return false;
2204 }
2205
2206 // deal with default values
2207 if(numPointsU < 1) {
2208 Msg::Error("Wrong number of control points along U for Bezier surface");
2209 return false;
2210 }
2211 int numPointsV = pointTags.size() / numPointsU;
2212 if(numPointsU * numPointsV != (int)pointTags.size()) {
2213 Msg::Error("Wrong number of control points for Bezier surface");
2214 return false;
2215 }
2216
2217 std::vector<TopoDS_Wire> wires;
2218 for(std::size_t i = 0; i < wireTags.size(); i++) {
2219 int wireTag = std::abs(wireTags[i]);
2220 if(!_tagWire.IsBound(wireTag)) {
2221 Msg::Error("Unknown OpenCASCADE curve loop with tag %d", wireTag);
2222 return false;
2223 }
2224 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTag));
2225 wires.push_back(wire);
2226 }
2227
2228 TopoDS_Face result;
2229 try {
2230 TColgp_Array2OfPnt pp(1, numPointsU, 1, numPointsV);
2231 for(int i = 1; i <= numPointsU; i++) {
2232 for(int j = 1; j <= numPointsV; j++) {
2233 int k = (j - 1) * numPointsU + (i - 1);
2234 if(!_tagVertex.IsBound(pointTags[k])) {
2235 Msg::Error("Unknown OpenCASCADE point with tag %d", pointTags[k]);
2236 return false;
2237 }
2238 TopoDS_Vertex vertex = TopoDS::Vertex(_tagVertex.Find(pointTags[k]));
2239 pp.SetValue(i, j, BRep_Tool::Pnt(vertex));
2240 }
2241 }
2242 Handle(Geom_BezierSurface) surf = new Geom_BezierSurface(pp);
2243 makeTrimmedSurface(surf, wires, wire3D, result);
2244 } catch(Standard_Failure &err) {
2245 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2246 return false;
2247 }
2248
2249 if(tag < 0) tag = getMaxTag(2) + 1;
2250 _bind(result, tag, true);
2251 return true;
2252 }
2253
addTrimmedSurface(int & tag,int surfaceTag,const std::vector<int> & wireTags,bool wire3D)2254 bool OCC_Internals::addTrimmedSurface(int &tag, int surfaceTag,
2255 const std::vector<int> &wireTags,
2256 bool wire3D)
2257 {
2258 if(tag >= 0 && _tagFace.IsBound(tag)) {
2259 Msg::Error("OpenCASCADE surface with tag %d already exists", tag);
2260 return false;
2261 }
2262
2263 if(!_tagFace.IsBound(surfaceTag)) {
2264 Msg::Error("Unknown OpenCASCADE surface with tag %d", surfaceTag);
2265 return false;
2266 }
2267 TopoDS_Face face = TopoDS::Face(_tagFace.Find(surfaceTag));
2268
2269 std::vector<TopoDS_Wire> wires;
2270 for(std::size_t i = 0; i < wireTags.size(); i++) {
2271 int wireTag = std::abs(wireTags[i]);
2272 if(!_tagWire.IsBound(wireTag)) {
2273 Msg::Error("Unknown OpenCASCADE curve loop with tag %d", wireTag);
2274 return false;
2275 }
2276 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTag));
2277 wires.push_back(wire);
2278 }
2279
2280 TopoDS_Face result;
2281 try {
2282 Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
2283 makeTrimmedSurface(surf, wires, wire3D, result);
2284 } catch(Standard_Failure &err) {
2285 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2286 return false;
2287 }
2288
2289 if(tag < 0) tag = getMaxTag(2) + 1;
2290 _bind(result, tag, true);
2291 return true;
2292 }
2293
addSurfaceLoop(int & tag,const std::vector<int> & surfaceTags,bool sewing)2294 bool OCC_Internals::addSurfaceLoop(int &tag,
2295 const std::vector<int> &surfaceTags,
2296 bool sewing)
2297 {
2298 if(tag >= 0 && _tagShell.IsBound(tag)) {
2299 Msg::Error("OpenCASCADE surface loop with tag %d already exists", tag);
2300 return false;
2301 }
2302
2303 if(sewing) {
2304 // this allows one to build a shell made of surfaces that share
2305 // geometrically identical (but topologically different) curves.
2306 TopoDS_Shape result;
2307 try {
2308 BRepBuilderAPI_Sewing s;
2309 for(std::size_t i = 0; i < surfaceTags.size(); i++) {
2310 if(!_tagFace.IsBound(surfaceTags[i])) {
2311 Msg::Error("Unknown OpenCASCADE surface with tag %d", surfaceTags[i]);
2312 return false;
2313 }
2314 TopoDS_Face face = TopoDS::Face(_tagFace.Find(surfaceTags[i]));
2315 s.Add(face);
2316 }
2317 s.Perform();
2318 result = s.SewedShape();
2319 } catch(Standard_Failure &err) {
2320 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2321 return false;
2322 }
2323 bool first = true;
2324 TopExp_Explorer exp0;
2325 for(exp0.Init(result, TopAbs_SHELL); exp0.More(); exp0.Next()) {
2326 TopoDS_Shell shell = TopoDS::Shell(exp0.Current());
2327 if(CTX::instance()->geom.occAutoFix) {
2328 // make sure faces in shell are oriented correctly
2329 ShapeFix_Shell fix(shell);
2330 fix.Perform();
2331 shell = fix.Shell();
2332 }
2333 int t = tag;
2334 if(first) { first = false; }
2335 else {
2336 t = getMaxTag(-2) + 1;
2337 Msg::Warning("Creating additional surface loop %d", t);
2338 }
2339 _bind(shell, t, true);
2340 return true;
2341 }
2342 }
2343
2344 try {
2345 BRep_Builder builder;
2346 BRepPrim_Builder b(builder);
2347 TopoDS_Shell shell;
2348 b.MakeShell(shell);
2349 for(std::size_t i = 0; i < surfaceTags.size(); i++) {
2350 if(!_tagFace.IsBound(surfaceTags[i])) {
2351 Msg::Error("Unknown OpenCASCADE surface with tag %d", surfaceTags[i]);
2352 return false;
2353 }
2354 TopoDS_Face face = TopoDS::Face(_tagFace.Find(surfaceTags[i]));
2355 b.AddShellFace(shell, face);
2356 }
2357 if(CTX::instance()->geom.occAutoFix) {
2358 // make sure faces in shell are oriented correctly
2359 ShapeFix_Shell fix(shell);
2360 fix.Perform();
2361 shell = fix.Shell();
2362 }
2363 _bind(shell, tag, true);
2364 return true;
2365 } catch(Standard_Failure &err) {
2366 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2367 return false;
2368 }
2369 }
2370
addVolume(int & tag,const std::vector<int> & shellTags)2371 bool OCC_Internals::addVolume(int &tag, const std::vector<int> &shellTags)
2372 {
2373 if(tag >= 0 && _tagSolid.IsBound(tag)) {
2374 Msg::Error("OpenCASCADE volume with tag %d already exists", tag);
2375 return false;
2376 }
2377
2378 TopoDS_Solid result;
2379 try {
2380 BRepBuilderAPI_MakeSolid s;
2381 for(std::size_t i = 0; i < shellTags.size(); i++) {
2382 if(!_tagShell.IsBound(shellTags[i])) {
2383 Msg::Error("Unknown OpenCASCADE surface loop with tag %d",
2384 shellTags[i]);
2385 return false;
2386 }
2387 TopoDS_Shell shell = TopoDS::Shell(_tagShell.Find(shellTags[i]));
2388 s.Add(shell);
2389 }
2390 result = s.Solid();
2391 if(CTX::instance()->geom.occAutoFix) {
2392 // make sure the volume is finite
2393 ShapeFix_Solid fix(result);
2394 fix.Perform();
2395 result = TopoDS::Solid(fix.Solid());
2396 }
2397 } catch(Standard_Failure &err) {
2398 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2399 return false;
2400 }
2401 if(tag < 0) tag = getMaxTag(3) + 1;
2402 _bind(result, tag, true);
2403 return true;
2404 }
2405
makeSphere(TopoDS_Solid & result,double xc,double yc,double zc,double radius,double angle1,double angle2,double angle3)2406 static bool makeSphere(TopoDS_Solid &result, double xc, double yc, double zc,
2407 double radius, double angle1, double angle2,
2408 double angle3)
2409 {
2410 if(radius <= 0) {
2411 Msg::Error("Sphere radius should be positive");
2412 return false;
2413 }
2414 if(angle3 <= 0 || angle3 > 2 * M_PI) {
2415 Msg::Error("Cannot build sphere with angle <= 0 or angle > 2*Pi");
2416 return false;
2417 }
2418 try {
2419 gp_Pnt p(xc, yc, zc);
2420 BRepPrimAPI_MakeSphere s(p, radius, angle1, angle2, angle3);
2421 s.Build();
2422 if(!s.IsDone()) {
2423 Msg::Error("Could not create sphere");
2424 return false;
2425 }
2426 result = TopoDS::Solid(s.Shape());
2427 } catch(Standard_Failure &err) {
2428 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2429 return false;
2430 }
2431 return true;
2432 }
2433
addSphere(int & tag,double xc,double yc,double zc,double radius,double angle1,double angle2,double angle3)2434 bool OCC_Internals::addSphere(int &tag, double xc, double yc, double zc,
2435 double radius, double angle1, double angle2,
2436 double angle3)
2437 {
2438 if(tag >= 0 && _tagSolid.IsBound(tag)) {
2439 Msg::Error("OpenCASCADE volume with tag %d already exists", tag);
2440 return false;
2441 }
2442 TopoDS_Solid result;
2443 if(!makeSphere(result, xc, yc, zc, radius, angle1, angle2, angle3))
2444 return false;
2445 if(tag < 0) tag = getMaxTag(3) + 1;
2446 _bind(result, tag, true);
2447 return true;
2448 }
2449
makeBox(TopoDS_Solid & result,double x,double y,double z,double dx,double dy,double dz)2450 static bool makeBox(TopoDS_Solid &result, double x, double y, double z,
2451 double dx, double dy, double dz)
2452 {
2453 if(!dx || !dy || !dz) {
2454 Msg::Error("Degenerate box");
2455 return false;
2456 }
2457 try {
2458 gp_Pnt P1(x, y, z);
2459 gp_Pnt P2(x + dx, y + dy, z + dz);
2460 BRepPrimAPI_MakeBox b(P1, P2);
2461 b.Build();
2462 if(!b.IsDone()) {
2463 Msg::Error("Could not create box");
2464 return false;
2465 }
2466 result = TopoDS::Solid(b.Shape());
2467 } catch(Standard_Failure &err) {
2468 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2469 return false;
2470 }
2471 return true;
2472 }
2473
addBox(int & tag,double x,double y,double z,double dx,double dy,double dz)2474 bool OCC_Internals::addBox(int &tag, double x, double y, double z, double dx,
2475 double dy, double dz)
2476 {
2477 if(tag >= 0 && _tagSolid.IsBound(tag)) {
2478 Msg::Error("OpenCASCADE volume with tag %d already exists", tag);
2479 return false;
2480 }
2481 TopoDS_Solid result;
2482 if(!makeBox(result, x, y, z, dx, dy, dz)) return false;
2483 if(tag < 0) tag = getMaxTag(3) + 1;
2484 _bind(result, tag, true);
2485 return true;
2486 }
2487
makeCylinder(TopoDS_Solid & result,double x,double y,double z,double dx,double dy,double dz,double r,double angle)2488 static bool makeCylinder(TopoDS_Solid &result, double x, double y, double z,
2489 double dx, double dy, double dz, double r,
2490 double angle)
2491 {
2492 const double H = sqrt(dx * dx + dy * dy + dz * dz);
2493 if(!H) {
2494 Msg::Error("Cannot build cylinder of zero height");
2495 return false;
2496 }
2497 if(angle <= 0 || angle > 2 * M_PI) {
2498 Msg::Error("Cannot build cylinder with angle <= 0 or angle > 2*Pi");
2499 return false;
2500 }
2501 try {
2502 gp_Pnt aP(x, y, z);
2503 gp_Vec aV(dx / H, dy / H, dz / H);
2504 gp_Ax2 anAxes(aP, aV);
2505 BRepPrimAPI_MakeCylinder c(anAxes, r, H, angle);
2506 c.Build();
2507 if(!c.IsDone()) {
2508 Msg::Error("Could not create cylinder");
2509 return false;
2510 }
2511 result = TopoDS::Solid(c.Shape());
2512 } catch(Standard_Failure &err) {
2513 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2514 return false;
2515 }
2516 return true;
2517 }
2518
addCylinder(int & tag,double x,double y,double z,double dx,double dy,double dz,double r,double angle)2519 bool OCC_Internals::addCylinder(int &tag, double x, double y, double z,
2520 double dx, double dy, double dz, double r,
2521 double angle)
2522 {
2523 if(tag >= 0 && _tagSolid.IsBound(tag)) {
2524 Msg::Error("OpenCASCADE volume with tag %d already exists", tag);
2525 return false;
2526 }
2527 TopoDS_Solid result;
2528 if(!makeCylinder(result, x, y, z, dx, dy, dz, r, angle)) return false;
2529 if(tag < 0) tag = getMaxTag(3) + 1;
2530 _bind(result, tag, true);
2531 return true;
2532 }
2533
makeTorus(TopoDS_Solid & result,double x,double y,double z,double r1,double r2,double angle)2534 static bool makeTorus(TopoDS_Solid &result, double x, double y, double z,
2535 double r1, double r2, double angle)
2536 {
2537 if(r1 <= 0 || r2 <= 0) {
2538 Msg::Error("Torus radii should be positive");
2539 return false;
2540 }
2541 try {
2542 gp_Pnt aP(x, y, z);
2543 gp_Vec aV(0, 0, 1);
2544 gp_Ax2 anAxes(aP, aV);
2545 BRepPrimAPI_MakeTorus t(anAxes, r1, r2, angle);
2546 t.Build();
2547 if(!t.IsDone()) {
2548 Msg::Error("Could not create torus");
2549 return false;
2550 }
2551 result = TopoDS::Solid(t.Shape());
2552 } catch(Standard_Failure &err) {
2553 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2554 return false;
2555 }
2556 return true;
2557 }
2558
addTorus(int & tag,double x,double y,double z,double r1,double r2,double angle)2559 bool OCC_Internals::addTorus(int &tag, double x, double y, double z, double r1,
2560 double r2, double angle)
2561 {
2562 if(tag >= 0 && _tagSolid.IsBound(tag)) {
2563 Msg::Error("OpenCASCADE volume with tag %d already exists", tag);
2564 return false;
2565 }
2566 TopoDS_Solid result;
2567 if(!makeTorus(result, x, y, z, r1, r2, angle)) return false;
2568 if(tag < 0) tag = getMaxTag(3) + 1;
2569 _bind(result, tag, true);
2570 return true;
2571 }
2572
makeCone(TopoDS_Solid & result,double x,double y,double z,double dx,double dy,double dz,double r1,double r2,double angle)2573 static bool makeCone(TopoDS_Solid &result, double x, double y, double z,
2574 double dx, double dy, double dz, double r1, double r2,
2575 double angle)
2576 {
2577 const double H = sqrt(dx * dx + dy * dy + dz * dz);
2578 if(!H) {
2579 Msg::Error("Cannot build cone of zero height");
2580 return false;
2581 }
2582 if(angle <= 0) {
2583 Msg::Error("Cone angle should be positive");
2584 return false;
2585 }
2586 try {
2587 gp_Pnt aP(x, y, z);
2588 gp_Vec aV(dx / H, dy / H, dz / H);
2589 gp_Ax2 anAxes(aP, aV);
2590 BRepPrimAPI_MakeCone c(anAxes, r1, r2, H, angle);
2591 c.Build();
2592 if(!c.IsDone()) {
2593 Msg::Error("Could not create cone");
2594 return false;
2595 }
2596 result = TopoDS::Solid(c.Shape());
2597 } catch(Standard_Failure &err) {
2598 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2599 return false;
2600 }
2601 return true;
2602 }
2603
addCone(int & tag,double x,double y,double z,double dx,double dy,double dz,double r1,double r2,double angle)2604 bool OCC_Internals::addCone(int &tag, double x, double y, double z, double dx,
2605 double dy, double dz, double r1, double r2,
2606 double angle)
2607 {
2608 if(tag >= 0 && _tagSolid.IsBound(tag)) {
2609 Msg::Error("OpenCASCADE volume with tag %d already exists", tag);
2610 return false;
2611 }
2612 TopoDS_Solid result;
2613 if(!makeCone(result, x, y, z, dx, dy, dz, r1, r2, angle)) return false;
2614 if(tag < 0) tag = getMaxTag(3) + 1;
2615 _bind(result, tag, true);
2616 return true;
2617 }
2618
makeWedge(TopoDS_Solid & result,double x,double y,double z,double dx,double dy,double dz,double ltx)2619 static bool makeWedge(TopoDS_Solid &result, double x, double y, double z,
2620 double dx, double dy, double dz, double ltx)
2621 {
2622 try {
2623 gp_Pnt aP(x, y, z);
2624 gp_Vec aV(0, 0, 1);
2625 gp_Ax2 anAxes(aP, aV);
2626 BRepPrimAPI_MakeWedge w(anAxes, dx, dy, dz, ltx);
2627 w.Build();
2628 if(!w.IsDone()) {
2629 Msg::Error("Could not create wedge");
2630 return false;
2631 }
2632 result = TopoDS::Solid(w.Shape());
2633 } catch(Standard_Failure &err) {
2634 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2635 return false;
2636 }
2637 return true;
2638 }
2639
addWedge(int & tag,double x,double y,double z,double dx,double dy,double dz,double ltx)2640 bool OCC_Internals::addWedge(int &tag, double x, double y, double z, double dx,
2641 double dy, double dz, double ltx)
2642 {
2643 if(tag >= 0 && _tagSolid.IsBound(tag)) {
2644 Msg::Error("OpenCASCADE volume with tag %d already exists", tag);
2645 return false;
2646 }
2647 TopoDS_Solid result;
2648 if(!makeWedge(result, x, y, z, dx, dy, dz, ltx)) return false;
2649 if(tag < 0) tag = getMaxTag(3) + 1;
2650 _bind(result, tag, true);
2651 return true;
2652 }
2653
addThruSections(int tag,const std::vector<int> & wireTags,bool makeSolid,bool makeRuled,std::vector<std::pair<int,int>> & outDimTags,int maxDegree)2654 bool OCC_Internals::addThruSections(
2655 int tag, const std::vector<int> &wireTags, bool makeSolid, bool makeRuled,
2656 std::vector<std::pair<int, int> > &outDimTags, int maxDegree)
2657 {
2658 int dim = makeSolid ? 3 : 2;
2659 if(tag >= 0 && _isBound(dim, tag)) {
2660 Msg::Error("OpenCASCADE entity of dimension %d with tag %d already exists",
2661 dim, tag);
2662 return false;
2663 }
2664 if(wireTags.size() < 2) {
2665 Msg::Error("ThruSections require at least 2 wires");
2666 return false;
2667 }
2668 TopoDS_Shape result;
2669 try {
2670 BRepOffsetAPI_ThruSections ts(makeSolid, makeRuled);
2671 // ts.SetContinuity(GeomAbs_C1);
2672 // Available choices:
2673 // GeomAbs_C0, GeomAbs_G1, GeomAbs_C1, GeomAbs_G2, GeomAbs_C2,
2674 // GeomAbs_C3, GeomAbs_CN
2675
2676 // ts.SetCriteriumWeight(1, 1, 1);
2677
2678 if(maxDegree > 0)
2679 ts.SetMaxDegree(maxDegree);
2680 else if(CTX::instance()->geom.occThruSectionsDegree > 0)
2681 ts.SetMaxDegree(CTX::instance()->geom.occThruSectionsDegree);
2682
2683 // ts.SetParType(Approx_ChordLength);
2684 // Available choices:
2685 // Approx_ChordLength, Approx_Centripetal, Approx_IsoParametric
2686
2687 // ts.SetSmoothing(Standard_True);
2688 for(std::size_t i = 0; i < wireTags.size(); i++) {
2689 if(!_tagWire.IsBound(wireTags[i])) {
2690 Msg::Error("Unknown OpenCASCADE wire or curve loop with tag %d",
2691 wireTags[i]);
2692 return false;
2693 }
2694 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTags[i]));
2695 if(makeSolid && !wire.Closed()) {
2696 Msg::Error("Making solid requires closed wires");
2697 return false;
2698 }
2699 ts.AddWire(wire);
2700 }
2701 ts.CheckCompatibility(Standard_False);
2702 ts.Build();
2703 if(!ts.IsDone()) {
2704 Msg::Error("Could not create ThruSection");
2705 return false;
2706 }
2707 result = ts.Shape();
2708 } catch(Standard_Failure &err) {
2709 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2710 return false;
2711 }
2712 _multiBind(result, tag, outDimTags, true, true);
2713 return true;
2714 }
2715
addThickSolid(int tag,int solidTag,const std::vector<int> & excludeFaceTags,double offset,std::vector<std::pair<int,int>> & outDimTags)2716 bool OCC_Internals::addThickSolid(int tag, int solidTag,
2717 const std::vector<int> &excludeFaceTags,
2718 double offset,
2719 std::vector<std::pair<int, int> > &outDimTags)
2720 {
2721 if(tag >= 0 && _isBound(3, tag)) {
2722 Msg::Error("OpenCASCADE volume with tag %d already exists", tag);
2723 return false;
2724 }
2725 if(!_isBound(3, solidTag)) {
2726 Msg::Error("Unknown OpenCASCADE volume with tag %d", solidTag);
2727 return false;
2728 }
2729 TopoDS_Shape result;
2730 try {
2731 TopoDS_Shape shape = _find(3, solidTag);
2732 TopTools_ListOfShape exclude;
2733 for(std::size_t i = 0; i < excludeFaceTags.size(); i++) {
2734 if(!_tagFace.IsBound(excludeFaceTags[i])) {
2735 Msg::Error("Unknown OpenCASCADE surface with tag %d",
2736 excludeFaceTags[i]);
2737 return false;
2738 }
2739 exclude.Append(_tagFace.Find(excludeFaceTags[i]));
2740 }
2741 #if OCC_VERSION_HEX > 0x070400
2742 BRepOffsetAPI_MakeThickSolid ts;
2743 ts.MakeThickSolidByJoin(shape, exclude, offset,
2744 CTX::instance()->geom.tolerance);
2745 #else
2746 BRepOffsetAPI_MakeThickSolid ts(shape, exclude, offset,
2747 CTX::instance()->geom.tolerance);
2748 ts.Build();
2749 #endif
2750 if(!ts.IsDone()) {
2751 Msg::Error("Could not build thick solid");
2752 return false;
2753 }
2754 result = ts.Shape();
2755 } catch(Standard_Failure &err) {
2756 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
2757 return false;
2758 }
2759 _multiBind(result, tag, outDimTags, true, true);
2760 return true;
2761 }
2762
_setExtrudedAttributes(const TopoDS_Compound & c,BRepSweep_Prism * p,BRepSweep_Revol * r,ExtrudeParams * e,double x,double y,double z,double dx,double dy,double dz,double ax,double ay,double az,double angle)2763 void OCC_Internals::_setExtrudedAttributes(
2764 const TopoDS_Compound &c, BRepSweep_Prism *p, BRepSweep_Revol *r,
2765 ExtrudeParams *e, double x, double y, double z, double dx, double dy,
2766 double dz, double ax, double ay, double az, double angle)
2767 {
2768 if(!p && !r) return;
2769
2770 bool extrude_attributes = (e ? true : false);
2771
2772 if(extrude_attributes && r && angle >= 2 * M_PI) {
2773 // OCC removes the origin edge from e.g. disks, which makes it impossible to
2774 // generate the 2D surface mesh by extrusion of the 1D edge mesh
2775 Msg::Warning("Extruded meshes by revolution only for angle < 2*Pi");
2776 extrude_attributes = false;
2777 }
2778
2779 TopExp_Explorer exp0;
2780
2781 for(exp0.Init(c, TopAbs_FACE); exp0.More(); exp0.Next()) {
2782 TopoDS_Face face = TopoDS::Face(exp0.Current());
2783 TopoDS_Shape bot = p ? p->FirstShape(face) : r->FirstShape(face);
2784 TopoDS_Shape top = p ? p->LastShape(face) : r->LastShape(face);
2785 if(extrude_attributes) {
2786 ExtrudeParams *ee = new ExtrudeParams(COPIED_ENTITY);
2787 ee->fill(p ? TRANSLATE : ROTATE, dx, dy, dz, ax, ay, az, x, y, z, angle);
2788 ee->mesh = e->mesh;
2789 _attributes->insert(new OCCAttributes(2, top, ee, 2, bot));
2790 }
2791 TopoDS_Shape vol = p ? p->Shape(face) : r->Shape(face);
2792 if(extrude_attributes) {
2793 ExtrudeParams *ee = new ExtrudeParams(EXTRUDED_ENTITY);
2794 ee->fill(p ? TRANSLATE : ROTATE, dx, dy, dz, ax, ay, az, x, y, z, angle);
2795 ee->mesh = e->mesh;
2796 _attributes->insert(new OCCAttributes(3, vol, ee, 2, bot));
2797 }
2798 }
2799
2800 for(exp0.Init(c, TopAbs_EDGE); exp0.More(); exp0.Next()) {
2801 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
2802 TopoDS_Shape bot = p ? p->FirstShape(edge) : r->FirstShape(edge);
2803 TopoDS_Shape top = p ? p->LastShape(edge) : r->LastShape(edge);
2804 if(extrude_attributes) {
2805 ExtrudeParams *ee = new ExtrudeParams(COPIED_ENTITY);
2806 ee->fill(p ? TRANSLATE : ROTATE, dx, dy, dz, ax, ay, az, x, y, z, angle);
2807 ee->mesh = e->mesh;
2808 _attributes->insert(new OCCAttributes(1, top, ee, 1, bot));
2809 }
2810 TopoDS_Shape sur = p ? p->Shape(edge) : r->Shape(edge);
2811 if(extrude_attributes) {
2812 ExtrudeParams *ee = new ExtrudeParams(EXTRUDED_ENTITY);
2813 ee->fill(p ? TRANSLATE : ROTATE, dx, dy, dz, ax, ay, az, x, y, z, angle);
2814 ee->mesh = e->mesh;
2815 _attributes->insert(new OCCAttributes(2, sur, ee, 1, bot));
2816 }
2817 }
2818
2819 for(exp0.Init(c, TopAbs_VERTEX); exp0.More(); exp0.Next()) {
2820 TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
2821 TopoDS_Shape bot = p ? p->FirstShape(vertex) : r->FirstShape(vertex);
2822 TopoDS_Shape top = p ? p->LastShape(vertex) : r->LastShape(vertex);
2823 TopoDS_Shape lin = p ? p->Shape(vertex) : r->Shape(vertex);
2824 if(extrude_attributes) {
2825 ExtrudeParams *ee = new ExtrudeParams(EXTRUDED_ENTITY);
2826 ee->fill(p ? TRANSLATE : ROTATE, dx, dy, dz, ax, ay, az, x, y, z, angle);
2827 ee->mesh = e->mesh;
2828 _attributes->insert(new OCCAttributes(1, lin, ee, 0, bot));
2829 }
2830 {
2831 double lc = _attributes->getMeshSize(0, bot);
2832 if(lc > 0 && lc < MAX_LC)
2833 _attributes->insert(new OCCAttributes(0, top, lc));
2834 }
2835 }
2836 }
2837
_getFuzzyTag(int dim,const TopoDS_Shape & s)2838 int OCC_Internals::_getFuzzyTag(int dim, const TopoDS_Shape &s)
2839 {
2840 if(_isBound(dim, s)) return _find(dim, s);
2841
2842 std::vector<TopoDS_Shape> candidates;
2843 _attributes->getSimilarShapes(dim, s, candidates);
2844
2845 int num = 0;
2846 for(std::size_t i = 0; i < candidates.size(); i++) {
2847 if(_isBound(dim, candidates[i])) { num++; }
2848 }
2849 Msg::Debug("Extruded mesh constraint fuzzy search: found %d candidates "
2850 "(dim=%d, %d bound)",
2851 (int)candidates.size(), dim, num);
2852 for(std::size_t i = 0; i < candidates.size(); i++) {
2853 if(_isBound(dim, candidates[i])) { return _find(dim, candidates[i]); }
2854 }
2855 return -1;
2856 }
2857
_copyExtrudedAttributes(TopoDS_Edge edge,GEdge * ge)2858 void OCC_Internals::_copyExtrudedAttributes(TopoDS_Edge edge, GEdge *ge)
2859 {
2860 int sourceDim = -1;
2861 TopoDS_Shape sourceShape;
2862 ExtrudeParams *e =
2863 _attributes->getExtrudeParams(1, edge, sourceDim, sourceShape);
2864 if(!e) return;
2865 if(e->geo.Mode == EXTRUDED_ENTITY) {
2866 e->geo.Source = _getFuzzyTag(0, sourceShape);
2867 }
2868 else if(e->geo.Mode == COPIED_ENTITY) {
2869 e->geo.Source = _getFuzzyTag(1, sourceShape);
2870 // detect degenerate extrusions or cycles
2871 ExtrudeParams *p = e;
2872 int recur = 0;
2873 while(++recur < CTX::instance()->mesh.maxRetries) {
2874 if(ge->tag() == p->geo.Source) {
2875 Msg::Info("Extrusion layer cycle detected for curve %d", ge->tag());
2876 e = nullptr;
2877 break;
2878 }
2879 GEdge *src = ge->model()->getEdgeByTag(p->geo.Source);
2880 if(src && src->meshAttributes.extrude &&
2881 src->meshAttributes.extrude->geo.Mode == COPIED_ENTITY) {
2882 p = src->meshAttributes.extrude;
2883 }
2884 else {
2885 break;
2886 }
2887 }
2888 }
2889 ge->meshAttributes.extrude = e;
2890 }
2891
_copyExtrudedAttributes(TopoDS_Face face,GFace * gf)2892 void OCC_Internals::_copyExtrudedAttributes(TopoDS_Face face, GFace *gf)
2893 {
2894 int sourceDim = -1;
2895 TopoDS_Shape sourceShape;
2896 ExtrudeParams *e =
2897 _attributes->getExtrudeParams(2, face, sourceDim, sourceShape);
2898 if(!e) return;
2899 if(e->geo.Mode == EXTRUDED_ENTITY) {
2900 e->geo.Source = _getFuzzyTag(1, sourceShape);
2901 }
2902 else if(e->geo.Mode == COPIED_ENTITY) {
2903 e->geo.Source = _getFuzzyTag(2, sourceShape);
2904 // detect degenerate extrusions or cycles
2905 ExtrudeParams *p = e;
2906 int recur = 0;
2907 while(++recur < CTX::instance()->mesh.maxRetries) {
2908 if(gf->tag() == p->geo.Source) {
2909 Msg::Info("Extrusion layer cycle detected for surface %d", gf->tag());
2910 e = nullptr;
2911 break;
2912 }
2913 GFace *src = gf->model()->getFaceByTag(p->geo.Source);
2914 if(src && src->meshAttributes.extrude &&
2915 src->meshAttributes.extrude->geo.Mode == COPIED_ENTITY) {
2916 p = src->meshAttributes.extrude;
2917 }
2918 else {
2919 break;
2920 }
2921 }
2922 }
2923 gf->meshAttributes.extrude = e;
2924 }
2925
_copyExtrudedAttributes(TopoDS_Solid solid,GRegion * gr)2926 void OCC_Internals::_copyExtrudedAttributes(TopoDS_Solid solid, GRegion *gr)
2927 {
2928 int sourceDim = -1;
2929 TopoDS_Shape sourceShape;
2930 ExtrudeParams *e =
2931 _attributes->getExtrudeParams(3, solid, sourceDim, sourceShape);
2932 if(!e) return;
2933 if(e->geo.Mode == EXTRUDED_ENTITY) {
2934 e->geo.Source = _getFuzzyTag(2, sourceShape);
2935 }
2936 gr->meshAttributes.extrude = e;
2937 }
2938
2939 template <class T>
getReturnedShapes(const TopoDS_Compound & c,T * sweep,std::vector<TopoDS_Shape> & top,std::vector<TopoDS_Shape> & body,std::vector<std::vector<TopoDS_Shape>> & lateral)2940 static int getReturnedShapes(const TopoDS_Compound &c, T *sweep,
2941 std::vector<TopoDS_Shape> &top,
2942 std::vector<TopoDS_Shape> &body,
2943 std::vector<std::vector<TopoDS_Shape> > &lateral)
2944 {
2945 TopExp_Explorer exp0, exp1;
2946 for(exp0.Init(c, TopAbs_FACE); exp0.More(); exp0.Next()) {
2947 TopoDS_Face face = TopoDS::Face(exp0.Current());
2948 top.push_back(sweep->LastShape(face));
2949 body.push_back(sweep->Shape(face));
2950 lateral.push_back(std::vector<TopoDS_Shape>());
2951 for(exp1.Init(face, TopAbs_EDGE); exp1.More(); exp1.Next()) {
2952 TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
2953 lateral.back().push_back(sweep->Shape(edge));
2954 }
2955 }
2956 if(top.size()) return 3;
2957 for(exp0.Init(c, TopAbs_EDGE); exp0.More(); exp0.Next()) {
2958 TopoDS_Edge edge = TopoDS::Edge(exp0.Current());
2959 top.push_back(sweep->LastShape(edge));
2960 body.push_back(sweep->Shape(edge));
2961 lateral.push_back(std::vector<TopoDS_Shape>());
2962 for(exp1.Init(edge, TopAbs_VERTEX); exp1.More(); exp1.Next()) {
2963 TopoDS_Vertex vertex = TopoDS::Vertex(exp1.Current());
2964 lateral.back().push_back(sweep->Shape(vertex));
2965 }
2966 }
2967 if(top.size()) return 2;
2968 for(exp0.Init(c, TopAbs_VERTEX); exp0.More(); exp0.Next()) {
2969 TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
2970 top.push_back(sweep->LastShape(vertex));
2971 body.push_back(sweep->Shape(vertex));
2972 }
2973 if(top.size()) return 1;
2974 return 0;
2975 }
2976
_extrudePerDim(int mode,int inDim,const std::vector<int> & inTags,double x,double y,double z,double dx,double dy,double dz,double ax,double ay,double az,double angle,int wireTag,std::vector<std::pair<int,int>> & outDimTags,ExtrudeParams * e,const std::string & trihedron)2977 bool OCC_Internals::_extrudePerDim(
2978 int mode, int inDim, const std::vector<int> &inTags, double x, double y,
2979 double z, double dx, double dy, double dz, double ax, double ay, double az,
2980 double angle, int wireTag, std::vector<std::pair<int, int> > &outDimTags,
2981 ExtrudeParams *e, const std::string &trihedron)
2982 {
2983 // build a single compound shape, so that we won't duplicate internal
2984 // boundaries
2985 BRep_Builder b;
2986 TopoDS_Compound c;
2987 b.MakeCompound(c);
2988 for(std::size_t i = 0; i < inTags.size(); i++) {
2989 if(!_isBound(inDim, inTags[i])) {
2990 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d",
2991 inDim, inTags[i]);
2992 return false;
2993 }
2994 TopoDS_Shape shape = _find(inDim, inTags[i]);
2995 b.Add(c, shape);
2996 }
2997 TopoDS_Shape result;
2998 std::vector<TopoDS_Shape> top, body;
2999 std::vector<std::vector<TopoDS_Shape> > lateral;
3000 int dim = -1;
3001 try {
3002 if(mode == 0) { // extrude
3003 BRepPrimAPI_MakePrism p(c, gp_Vec(dx, dy, dz), Standard_False);
3004 p.Build();
3005 if(!p.IsDone()) {
3006 Msg::Error("Could not extrude");
3007 return false;
3008 }
3009 result = p.Shape();
3010 const BRepSweep_Prism &prism(p.Prism());
3011 _setExtrudedAttributes(c, (BRepSweep_Prism *)&prism, nullptr, e, 0., 0.,
3012 0., dx, dy, dz, 0., 0., 0., 0.);
3013 dim = getReturnedShapes(c, (BRepSweep_Prism *)&prism, top, body, lateral);
3014 }
3015 else if(mode == 1) { // revolve
3016 gp_Ax1 axisOfRevolution(gp_Pnt(x, y, z), gp_Dir(ax, ay, az));
3017 BRepPrimAPI_MakeRevol r(c, axisOfRevolution, angle, Standard_False);
3018 r.Build();
3019 if(!r.IsDone()) {
3020 Msg::Error("Could not revolve");
3021 return false;
3022 }
3023 result = r.Shape();
3024 const BRepSweep_Revol &revol(r.Revol());
3025 _setExtrudedAttributes(c, nullptr, (BRepSweep_Revol *)&revol, e, x, y, z,
3026 0., 0., 0., ax, ay, az, angle);
3027 dim = getReturnedShapes(c, (BRepSweep_Revol *)&revol, top, body, lateral);
3028 }
3029 else if(mode == 2) { // pipe
3030 if(!_tagWire.IsBound(wireTag)) {
3031 Msg::Error("Unknown OpenCASCADE wire with tag %d", wireTag);
3032 return false;
3033 }
3034 TopoDS_Wire wire = TopoDS::Wire(_tagWire.Find(wireTag));
3035 // DiscreteTrihedron seems the most robust; CorrectedFrenet e.g. fails on
3036 // very simple cases with straight extrusions.
3037 GeomFill_Trihedron mode = GeomFill_IsDiscreteTrihedron;
3038 if(trihedron == "" || trihedron == "DiscreteTrihedron")
3039 mode = GeomFill_IsDiscreteTrihedron;
3040 else if(trihedron == "CorrectedFrenet")
3041 mode = GeomFill_IsCorrectedFrenet;
3042 else if(trihedron == "Fixed")
3043 mode = GeomFill_IsFixed;
3044 else if(trihedron == "Frenet")
3045 mode = GeomFill_IsFrenet;
3046 else if(trihedron == "ConstantNormal")
3047 mode = GeomFill_IsConstantNormal;
3048 else if(trihedron == "Darboux")
3049 mode = GeomFill_IsDarboux;
3050 else if(trihedron == "GuideAC")
3051 mode = GeomFill_IsGuideAC;
3052 else if(trihedron == "GuidePlan")
3053 mode = GeomFill_IsGuidePlan;
3054 else if(trihedron == "GuideACWithContact")
3055 mode = GeomFill_IsGuideACWithContact;
3056 else if(trihedron == "GuidePlanWithContact")
3057 mode = GeomFill_IsGuidePlanWithContact;
3058 else
3059 Msg::Warning(
3060 "Unknown trihedron mode for pipe: using 'DiscreteTrihedron'");
3061 BRepOffsetAPI_MakePipe p(wire, c, mode);
3062 p.Build();
3063 if(!p.IsDone()) {
3064 Msg::Error("Could not create pipe");
3065 return false;
3066 }
3067 result = p.Shape();
3068 // const BRepFill_Pipe &pipe(p.Pipe());
3069 if(e)
3070 Msg::Warning(
3071 "Structured meshes not yet available with OpenCASCADE pipe");
3072 // Check if
3073 // pipe.FirstShape() gives us "bottom"
3074 // pipe.LastShape() gives us "top"
3075 // pipe.Shape() gives us "body"
3076 // using pipe.Spine(), pipe.{Face,Edge}(spine, c) gives us the lateral
3077 // entities
3078 // dim = getReturnedShapesForPipe(c, pipe, top, body, lateral);
3079 }
3080 } catch(Standard_Failure &err) {
3081 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
3082 return false;
3083 }
3084
3085 _multiBind(result, -1, outDimTags, true, true);
3086
3087 // return entities in the same order as the built-in kernel extrusion
3088 if(dim >= 1 && dim <= 3 && top.size() == inTags.size() &&
3089 top.size() == body.size()) {
3090 outDimTags.clear();
3091 for(std::size_t i = 0; i < top.size(); i++) {
3092 if(_isBound(dim - 1, top[i]))
3093 outDimTags.push_back(
3094 std::make_pair(dim - 1, _find(dim - 1, top[i])));
3095 if(_isBound(dim, body[i]))
3096 outDimTags.push_back(std::make_pair(dim, _find(dim, body[i])));
3097 if(CTX::instance()->geom.extrudeReturnLateral &&
3098 top.size() == lateral.size()) {
3099 for(std::size_t j = 0; j < lateral[i].size(); j++) {
3100 if(_isBound(dim - 1, lateral[i][j]))
3101 outDimTags.push_back(
3102 std::make_pair(dim - 1, _find(dim - 1, lateral[i][j])));
3103 }
3104 }
3105 }
3106 }
3107 return true;
3108 }
3109
_extrude(int mode,const std::vector<std::pair<int,int>> & inDimTags,double x,double y,double z,double dx,double dy,double dz,double ax,double ay,double az,double angle,int wireTag,std::vector<std::pair<int,int>> & outDimTags,ExtrudeParams * e,const std::string & trihedron)3110 bool OCC_Internals::_extrude(int mode,
3111 const std::vector<std::pair<int, int> > &inDimTags,
3112 double x, double y, double z, double dx, double dy,
3113 double dz, double ax, double ay, double az,
3114 double angle, int wireTag,
3115 std::vector<std::pair<int, int> > &outDimTags,
3116 ExtrudeParams *e, const std::string &trihedron)
3117 {
3118 std::vector<int> inTags[4];
3119 for(std::size_t i = 0; i < inDimTags.size(); i++) {
3120 int dim = inDimTags[i].first;
3121 if(dim < 0 || dim > 3) {
3122 Msg::Error("Wrong input dimension in extrusion");
3123 return false;
3124 }
3125 inTags[dim].push_back(inDimTags[i].second);
3126 }
3127 for(int dim = 0; dim < 4; dim++) {
3128 if(!inTags[dim].empty()) {
3129 std::vector<std::pair<int, int> > out;
3130 if(_extrudePerDim(mode, dim, inTags[dim], x, y, z, dx, dy, dz, ax, ay, az,
3131 angle, wireTag, out, e, trihedron)) {
3132 outDimTags.insert(outDimTags.end(), out.begin(), out.end());
3133 }
3134 }
3135 }
3136 return true;
3137 }
3138
extrude(const std::vector<std::pair<int,int>> & inDimTags,double dx,double dy,double dz,std::vector<std::pair<int,int>> & outDimTags,ExtrudeParams * e)3139 bool OCC_Internals::extrude(const std::vector<std::pair<int, int> > &inDimTags,
3140 double dx, double dy, double dz,
3141 std::vector<std::pair<int, int> > &outDimTags,
3142 ExtrudeParams *e)
3143 {
3144 return _extrude(0, inDimTags, 0., 0., 0., dx, dy, dz, 0., 0., 0., 0., 0,
3145 outDimTags, e);
3146 }
3147
revolve(const std::vector<std::pair<int,int>> & inDimTags,double x,double y,double z,double ax,double ay,double az,double angle,std::vector<std::pair<int,int>> & outDimTags,ExtrudeParams * e)3148 bool OCC_Internals::revolve(const std::vector<std::pair<int, int> > &inDimTags,
3149 double x, double y, double z, double ax, double ay,
3150 double az, double angle,
3151 std::vector<std::pair<int, int> > &outDimTags,
3152 ExtrudeParams *e)
3153 {
3154 return _extrude(1, inDimTags, x, y, z, 0., 0., 0., ax, ay, az, angle, 0,
3155 outDimTags, e);
3156 }
3157
addPipe(const std::vector<std::pair<int,int>> & inDimTags,int wireTag,std::vector<std::pair<int,int>> & outDimTags,const std::string & trihedron)3158 bool OCC_Internals::addPipe(const std::vector<std::pair<int, int> > &inDimTags,
3159 int wireTag,
3160 std::vector<std::pair<int, int> > &outDimTags,
3161 const std::string &trihedron)
3162 {
3163 return _extrude(2, inDimTags, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., wireTag,
3164 outDimTags, nullptr, trihedron);
3165 }
3166
_fillet(int mode,const std::vector<int> & volumeTags,const std::vector<int> & curveTags,const std::vector<int> & surfaceTags,const std::vector<double> & param,std::vector<std::pair<int,int>> & outDimTags,bool removeVolume)3167 bool OCC_Internals::_fillet(int mode, const std::vector<int> &volumeTags,
3168 const std::vector<int> &curveTags,
3169 const std::vector<int> &surfaceTags,
3170 const std::vector<double> ¶m,
3171 std::vector<std::pair<int, int> > &outDimTags,
3172 bool removeVolume)
3173 {
3174 std::vector<TopoDS_Edge> edges;
3175 for(std::size_t i = 0; i < curveTags.size(); i++) {
3176 if(!_tagEdge.IsBound(curveTags[i])) {
3177 Msg::Error("Unknown OpenCASCADE curve with tag %d", curveTags[i]);
3178 return false;
3179 }
3180 edges.push_back(TopoDS::Edge(_tagEdge.Find(curveTags[i])));
3181 }
3182
3183 std::vector<TopoDS_Face> faces;
3184 for(std::size_t i = 0; i < surfaceTags.size(); i++) {
3185 if(!_tagFace.IsBound(surfaceTags[i])) {
3186 Msg::Error("Unknown OpenCASCADE surface with tag %d", surfaceTags[i]);
3187 return false;
3188 }
3189 faces.push_back(TopoDS::Face(_tagFace.Find(surfaceTags[i])));
3190 }
3191 if(mode && edges.size() != faces.size()) {
3192 Msg::Error("Different number of curves and surfaces for chamfer");
3193 return false;
3194 }
3195
3196 // build a single compound shape
3197 BRep_Builder b;
3198 TopoDS_Compound c;
3199 b.MakeCompound(c);
3200 for(std::size_t i = 0; i < volumeTags.size(); i++) {
3201 if(!_isBound(3, volumeTags[i])) {
3202 Msg::Error("Unknown OpenCASCADE volume with tag %d", volumeTags[i]);
3203 return false;
3204 }
3205 TopoDS_Shape shape = _find(3, volumeTags[i]);
3206 if(removeVolume) _unbind(shape, 3, volumeTags[i], true);
3207 if(CTX::instance()->geom.occAutoFix) {
3208 // make sure the volume is finite
3209 ShapeFix_Solid fix(TopoDS::Solid(shape));
3210 fix.Perform();
3211 shape = fix.Solid();
3212 }
3213 b.Add(c, shape);
3214 }
3215 TopoDS_Shape result;
3216 try {
3217 if(mode == 0) { // fillet
3218 BRepFilletAPI_MakeFillet f(c);
3219 for(std::size_t i = 0; i < edges.size(); i++) {
3220 if(param.size() == 1)
3221 f.Add(param[0], edges[i]);
3222 else if(param.size() == edges.size())
3223 f.Add(param[i], edges[i]);
3224 else if(param.size() == 2 * edges.size())
3225 f.Add(param[2 * i], param[2 * i + 1], edges[i]);
3226 }
3227 f.Build();
3228 if(!f.IsDone()) {
3229 Msg::Error("Could not compute fillet");
3230 return false;
3231 }
3232 result = f.Shape();
3233 }
3234 else { // chamfer
3235 BRepFilletAPI_MakeChamfer f(c);
3236 for(std::size_t i = 0; i < edges.size(); i++) {
3237 if(param.size() == 1)
3238 f.Add(param[0], param[0], edges[i], faces[i]);
3239 else if(param.size() == edges.size())
3240 f.Add(param[i], param[i], edges[i], faces[i]);
3241 else if(param.size() == 2 * edges.size())
3242 f.Add(param[2 * i], param[2 * i + 1], edges[i], faces[i]);
3243 }
3244 f.Build();
3245 if(!f.IsDone()) {
3246 Msg::Error("Could not compute chamfer");
3247 return false;
3248 }
3249 result = f.Shape();
3250 }
3251 } catch(Standard_Failure &err) {
3252 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
3253 return false;
3254 }
3255
3256 if(result.IsNull()) {
3257 Msg::Error("%s produces empty shape", mode ? "Chamfer" : "Fillet");
3258 return false;
3259 }
3260
3261 // TODO: if removeVolume and CTX::instance()->geom.occBooleanPreserveNumbering
3262 // are set we could use Generated(), Modified() and IsDeleted() in a similar
3263 // way as what we do for boolean operation, in order to try to preserve tags
3264
3265 _multiBind(result, -1, outDimTags, true, true);
3266 return true;
3267 }
3268
fillet(const std::vector<int> & volumeTags,const std::vector<int> & curveTags,const std::vector<double> & radii,std::vector<std::pair<int,int>> & outDimTags,bool removeVolume)3269 bool OCC_Internals::fillet(const std::vector<int> &volumeTags,
3270 const std::vector<int> &curveTags,
3271 const std::vector<double> &radii,
3272 std::vector<std::pair<int, int> > &outDimTags,
3273 bool removeVolume)
3274 {
3275 std::vector<int> dummy;
3276 return _fillet(0, volumeTags, curveTags, dummy, radii, outDimTags,
3277 removeVolume);
3278 }
3279
chamfer(const std::vector<int> & volumeTags,const std::vector<int> & curveTags,const std::vector<int> & surfaceTags,const std::vector<double> & distances,std::vector<std::pair<int,int>> & outDimTags,bool removeVolume)3280 bool OCC_Internals::chamfer(const std::vector<int> &volumeTags,
3281 const std::vector<int> &curveTags,
3282 const std::vector<int> &surfaceTags,
3283 const std::vector<double> &distances,
3284 std::vector<std::pair<int, int> > &outDimTags,
3285 bool removeVolume)
3286 {
3287 return _fillet(1, volumeTags, curveTags, surfaceTags, distances, outDimTags,
3288 removeVolume);
3289 }
3290
_filterTags(std::vector<std::pair<int,int>> & outDimTags,int minDim)3291 static void _filterTags(std::vector<std::pair<int, int> > &outDimTags,
3292 int minDim)
3293 {
3294 std::vector<std::pair<int, int> > tmp(outDimTags);
3295 outDimTags.clear();
3296 for(std::size_t i = 0; i < tmp.size(); i++) {
3297 if(tmp[i].first >= minDim) outDimTags.push_back(tmp[i]);
3298 }
3299 }
3300
booleanOperator(int tag,BooleanOperator op,const std::vector<std::pair<int,int>> & objectDimTags,const std::vector<std::pair<int,int>> & toolDimTags,std::vector<std::pair<int,int>> & outDimTags,std::vector<std::vector<std::pair<int,int>>> & outDimTagsMap,bool removeObject,bool removeTool)3301 bool OCC_Internals::booleanOperator(
3302 int tag, BooleanOperator op,
3303 const std::vector<std::pair<int, int> > &objectDimTags,
3304 const std::vector<std::pair<int, int> > &toolDimTags,
3305 std::vector<std::pair<int, int> > &outDimTags,
3306 std::vector<std::vector<std::pair<int, int> > > &outDimTagsMap,
3307 bool removeObject, bool removeTool)
3308 {
3309 double tolerance = CTX::instance()->geom.toleranceBoolean;
3310 bool parallel = CTX::instance()->geom.occParallel;
3311 bool preserveNumbering = CTX::instance()->geom.occBooleanPreserveNumbering;
3312
3313 if(objectDimTags.empty()) return true;
3314
3315 if(tag >= 0 && _isBound(objectDimTags[0].first, tag)) {
3316 Msg::Error("OpenCASCADE entity with tag %d already exists", tag);
3317 return false;
3318 }
3319
3320 int minDim = 3;
3321 TopTools_ListOfShape objectShapes, toolShapes;
3322 for(std::size_t i = 0; i < objectDimTags.size(); i++) {
3323 int dim = objectDimTags[i].first;
3324 int t = objectDimTags[i].second;
3325 if(!_isBound(dim, t)) {
3326 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
3327 t);
3328 return false;
3329 }
3330 else {
3331 TopoDS_Shape object = _find(dim, t);
3332 objectShapes.Append(object);
3333 }
3334 minDim = std::min(minDim, dim);
3335 }
3336 for(std::size_t i = 0; i < toolDimTags.size(); i++) {
3337 int dim = toolDimTags[i].first;
3338 int t = toolDimTags[i].second;
3339 if(!_isBound(dim, t)) {
3340 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
3341 t);
3342 return false;
3343 }
3344 else {
3345 TopoDS_Shape tool = _find(dim, t);
3346 toolShapes.Append(tool);
3347 }
3348 minDim = std::min(minDim, dim);
3349 }
3350
3351 TopoDS_Shape result;
3352 std::vector<TopoDS_Shape> mapOriginal;
3353 std::vector<TopTools_ListOfShape> mapModified, mapGenerated;
3354 std::vector<bool> mapDeleted;
3355 try {
3356 switch(op) {
3357 case OCC_Internals::Union: {
3358 BRepAlgoAPI_Fuse fuse;
3359 fuse.SetRunParallel(parallel);
3360 fuse.SetArguments(objectShapes);
3361 fuse.SetTools(toolShapes);
3362 if(tolerance > 0.0) fuse.SetFuzzyValue(tolerance);
3363
3364 // TODO: add gluing option to speed-up operations when no "real"
3365 // intersections are present
3366 // * default:
3367 // fuse.SetGlue(BOPAlgo_GlueOff);
3368 // * speed up if no real intersection but partial or full overlapping
3369 // faces/edges:
3370 // fuse.SetGlue(BOPAlgo_Shift);
3371 // * speed up if no real intersection and no partial overlaps:
3372 // fuse.SetGlue(BOPAlgo_GlueFull);
3373
3374 // TODO: add option to prevent reuse of existing shapes:
3375 // fuse.SetNonDestructive(true);
3376
3377 fuse.Build();
3378 if(!fuse.IsDone()) {
3379 Msg::Error("Fuse operation cannot be performed");
3380 return false;
3381 }
3382 if(CTX::instance()->geom.occUnionUnify) {
3383 // try to unify faces and edges of the shape (remove internal seams)
3384 // which lie on the same geometry
3385 #if OCC_VERSION_HEX < 0x070400
3386 result = fuse.Shape();
3387 ShapeUpgrade_UnifySameDomain unify(result);
3388 unify.Build();
3389 result = unify.Shape();
3390 #else
3391 // better, as it preserves the history; TODO: maybe we should also make
3392 // this available for the other boolean operations
3393 fuse.SimplifyResult();
3394 result = fuse.Shape();
3395 #endif
3396 }
3397 else {
3398 result = fuse.Shape();
3399 }
3400 TopTools_ListIteratorOfListOfShape it(objectShapes);
3401 for(; it.More(); it.Next()) {
3402 mapOriginal.push_back(it.Value());
3403 mapModified.push_back(fuse.Modified(it.Value()));
3404 mapDeleted.push_back(fuse.IsDeleted(it.Value()));
3405 mapGenerated.push_back(fuse.Generated(it.Value()));
3406 }
3407 TopTools_ListIteratorOfListOfShape it2(toolShapes);
3408 for(; it2.More(); it2.Next()) {
3409 mapOriginal.push_back(it2.Value());
3410 mapModified.push_back(fuse.Modified(it2.Value()));
3411 mapDeleted.push_back(fuse.IsDeleted(it2.Value()));
3412 mapGenerated.push_back(fuse.Generated(it2.Value()));
3413 }
3414 } break;
3415 case OCC_Internals::Intersection: {
3416 BRepAlgoAPI_Common common;
3417 common.SetRunParallel(parallel);
3418 common.SetArguments(objectShapes);
3419 common.SetTools(toolShapes);
3420 if(tolerance > 0.0) common.SetFuzzyValue(tolerance);
3421 common.Build();
3422 if(!common.IsDone()) {
3423 Msg::Error("Intersection operation cannot be performed");
3424 return false;
3425 }
3426 result = common.Shape();
3427 TopTools_ListIteratorOfListOfShape it(objectShapes);
3428 for(; it.More(); it.Next()) {
3429 mapOriginal.push_back(it.Value());
3430 mapModified.push_back(common.Modified(it.Value()));
3431 mapDeleted.push_back(common.IsDeleted(it.Value()));
3432 mapGenerated.push_back(common.Generated(it.Value()));
3433 }
3434 TopTools_ListIteratorOfListOfShape it2(toolShapes);
3435 for(; it2.More(); it2.Next()) {
3436 mapOriginal.push_back(it2.Value());
3437 mapModified.push_back(common.Modified(it2.Value()));
3438 mapDeleted.push_back(common.IsDeleted(it2.Value()));
3439 mapGenerated.push_back(common.Generated(it2.Value()));
3440 }
3441 } break;
3442
3443 case OCC_Internals::Difference: {
3444 BRepAlgoAPI_Cut cut;
3445 cut.SetRunParallel(parallel);
3446 cut.SetArguments(objectShapes);
3447 cut.SetTools(toolShapes);
3448 if(tolerance > 0.0) cut.SetFuzzyValue(tolerance);
3449 cut.Build();
3450 if(!cut.IsDone()) {
3451 Msg::Error("Intersection operation cannot be performed");
3452 return false;
3453 }
3454 result = cut.Shape();
3455 TopTools_ListIteratorOfListOfShape it(objectShapes);
3456 for(; it.More(); it.Next()) {
3457 mapOriginal.push_back(it.Value());
3458 mapModified.push_back(cut.Modified(it.Value()));
3459 mapDeleted.push_back(cut.IsDeleted(it.Value()));
3460 mapGenerated.push_back(cut.Generated(it.Value()));
3461 }
3462 TopTools_ListIteratorOfListOfShape it2(toolShapes);
3463 for(; it2.More(); it2.Next()) {
3464 mapOriginal.push_back(it2.Value());
3465 mapModified.push_back(cut.Modified(it2.Value()));
3466 mapDeleted.push_back(cut.IsDeleted(it2.Value()));
3467 mapGenerated.push_back(cut.Generated(it2.Value()));
3468 }
3469 } break;
3470
3471 case OCC_Internals::Fragments:
3472 default: {
3473 BRepAlgoAPI_BuilderAlgo fragments;
3474 fragments.SetRunParallel(parallel);
3475 objectShapes.Append(toolShapes);
3476 toolShapes.Clear();
3477 fragments.SetArguments(objectShapes);
3478 if(tolerance > 0.0) fragments.SetFuzzyValue(tolerance);
3479 fragments.Build();
3480 if(!fragments.IsDone()) {
3481 Msg::Error("Boolean fragments failed");
3482 return false;
3483 }
3484 result = fragments.Shape();
3485 TopTools_ListIteratorOfListOfShape it(objectShapes);
3486 for(; it.More(); it.Next()) {
3487 mapOriginal.push_back(it.Value());
3488 mapModified.push_back(fragments.Modified(it.Value()));
3489 mapDeleted.push_back(fragments.IsDeleted(it.Value()));
3490 mapGenerated.push_back(fragments.Generated(it.Value()));
3491 }
3492 } break;
3493 }
3494 } catch(Standard_Failure &err) {
3495 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
3496 return false;
3497 }
3498
3499 std::vector<std::pair<int, int> > inDimTags;
3500 inDimTags.insert(inDimTags.end(), objectDimTags.begin(), objectDimTags.end());
3501 inDimTags.insert(inDimTags.end(), toolDimTags.begin(), toolDimTags.end());
3502 std::size_t numObjects = objectDimTags.size();
3503
3504 if(tag >= 0 || !preserveNumbering) {
3505 // if we specify the tag explicitly, or if we don't care about preserving
3506 // the numering, just go ahead and bind the resulting shape (and sub-shapes)
3507 for(std::size_t i = 0; i < inDimTags.size(); i++) {
3508 bool remove = (i < numObjects) ? removeObject : removeTool;
3509 if(remove) {
3510 int d = inDimTags[i].first;
3511 int t = inDimTags[i].second;
3512 if(_isBound(d, t)) _unbind(_find(d, t), d, t, true);
3513 }
3514 }
3515 _multiBind(result, tag, outDimTags, (tag >= 0) ? true : false, true,
3516 (tag >= 0) ? false : true);
3517 _filterTags(outDimTags, minDim);
3518 }
3519 else {
3520 // otherwise, try to preserve the numbering of the input shapes that did not
3521 // change, or that were replaced by a single shape. Note that to preserve
3522 // the numbering of smaller dimension entities (on boundaries) they should
3523 // appear *before* higher dimensional entities in the object/tool lists.
3524 _toPreserve.clear();
3525 for(std::size_t i = 0; i < inDimTags.size(); i++) {
3526 int dim = inDimTags[i].first;
3527 int tag = inDimTags[i].second;
3528 bool remove = (i < numObjects) ? removeObject : removeTool;
3529 if(mapDeleted[i]) { // deleted
3530 if(remove) _unbind(mapOriginal[i], dim, tag, true);
3531 Msg::Debug("BOOL (%d,%d) deleted", dim, tag);
3532 }
3533 else if(mapModified[i].Extent() == 0) { // not modified
3534 auto ins = _toPreserve.insert(std::make_pair(dim, tag));
3535 if(ins.second) // it's not yet in outDimTags
3536 outDimTags.push_back(std::make_pair(dim, tag));
3537 Msg::Debug("BOOL (%d,%d) not modified", dim, tag);
3538 }
3539 else if(mapModified[i].Extent() == 1) { // replaced by single one
3540 if(remove) {
3541 _unbind(mapOriginal[i], dim, tag, true);
3542 _bind(mapModified[i].First(), dim, tag, false); // not recursive!
3543 int t = _find(dim, mapModified[i].First());
3544 if(tag != t)
3545 Msg::Info("Could not preserve tag of %dD object %d (->%d)", dim,
3546 tag, t);
3547 auto ins = _toPreserve.insert(std::make_pair(dim, t));
3548 if(ins.second) // it's not yet in outDimTags
3549 outDimTags.push_back(std::make_pair(dim, t));
3550 }
3551 Msg::Debug("BOOL (%d,%d) replaced by 1", dim, tag);
3552 }
3553 else {
3554 if(remove) _unbind(mapOriginal[i], dim, tag, true);
3555 Msg::Debug("BOOL (%d,%d) other", dim, tag);
3556 }
3557 }
3558 for(int d = -2; d <= 3; d++) _recomputeMaxTag(d);
3559 // bind all remaining entities and add the new ones to the returned list
3560 _multiBind(result, -1, outDimTags, false, true, true);
3561 _filterTags(outDimTags, minDim);
3562 _toPreserve.clear();
3563 }
3564
3565 // return input/output correspondance maps
3566 for(std::size_t i = 0; i < inDimTags.size(); i++) {
3567 int dim = inDimTags[i].first;
3568 int tag = inDimTags[i].second;
3569 std::pair<int, int> dimTag(dim, tag);
3570 std::vector<std::pair<int, int> > dimTags;
3571 if(mapDeleted[i]) { // deleted
3572 }
3573 else if(mapModified[i].Extent() == 0) { // not modified
3574 if(_isBound(dim, tag)) dimTags.push_back(dimTag);
3575 }
3576 else {
3577 TopTools_ListIteratorOfListOfShape it(mapModified[i]);
3578 for(; it.More(); it.Next()) {
3579 if(_isBound(dim, it.Value())) {
3580 int t = _find(dim, it.Value());
3581 dimTags.push_back(std::make_pair(dim, t));
3582 }
3583 }
3584 TopTools_ListIteratorOfListOfShape it2(mapGenerated[i]);
3585 for(; it2.More(); it2.Next()) {
3586 if(_isBound(dim, it2.Value())) {
3587 int t = _find(dim, it2.Value());
3588 dimTags.push_back(std::make_pair(dim, t));
3589 }
3590 }
3591 }
3592 std::ostringstream sstream;
3593 sstream << "BOOL in (" << dim << "," << tag << ") -> out";
3594 for(std::size_t j = 0; j < dimTags.size(); j++)
3595 sstream << " (" << dimTags[j].first << "," << dimTags[j].second << ")";
3596 Msg::Debug("%s", sstream.str().c_str());
3597 outDimTagsMap.push_back(dimTags);
3598 }
3599
3600 return true;
3601 }
3602
booleanUnion(int tag,const std::vector<std::pair<int,int>> & objectDimTags,const std::vector<std::pair<int,int>> & toolDimTags,std::vector<std::pair<int,int>> & outDimTags,std::vector<std::vector<std::pair<int,int>>> & outDimTagsMap,bool removeObject,bool removeTool)3603 bool OCC_Internals::booleanUnion(
3604 int tag, const std::vector<std::pair<int, int> > &objectDimTags,
3605 const std::vector<std::pair<int, int> > &toolDimTags,
3606 std::vector<std::pair<int, int> > &outDimTags,
3607 std::vector<std::vector<std::pair<int, int> > > &outDimTagsMap,
3608 bool removeObject, bool removeTool)
3609 {
3610 return booleanOperator(tag, OCC_Internals::Union, objectDimTags, toolDimTags,
3611 outDimTags, outDimTagsMap, removeObject, removeTool);
3612 }
3613
booleanIntersection(int tag,const std::vector<std::pair<int,int>> & objectDimTags,const std::vector<std::pair<int,int>> & toolDimTags,std::vector<std::pair<int,int>> & outDimTags,std::vector<std::vector<std::pair<int,int>>> & outDimTagsMap,bool removeObject,bool removeTool)3614 bool OCC_Internals::booleanIntersection(
3615 int tag, const std::vector<std::pair<int, int> > &objectDimTags,
3616 const std::vector<std::pair<int, int> > &toolDimTags,
3617 std::vector<std::pair<int, int> > &outDimTags,
3618 std::vector<std::vector<std::pair<int, int> > > &outDimTagsMap,
3619 bool removeObject, bool removeTool)
3620 {
3621 return booleanOperator(tag, OCC_Internals::Intersection, objectDimTags,
3622 toolDimTags, outDimTags, outDimTagsMap, removeObject,
3623 removeTool);
3624 }
3625
booleanDifference(int tag,const std::vector<std::pair<int,int>> & objectDimTags,const std::vector<std::pair<int,int>> & toolDimTags,std::vector<std::pair<int,int>> & outDimTags,std::vector<std::vector<std::pair<int,int>>> & outDimTagsMap,bool removeObject,bool removeTool)3626 bool OCC_Internals::booleanDifference(
3627 int tag, const std::vector<std::pair<int, int> > &objectDimTags,
3628 const std::vector<std::pair<int, int> > &toolDimTags,
3629 std::vector<std::pair<int, int> > &outDimTags,
3630 std::vector<std::vector<std::pair<int, int> > > &outDimTagsMap,
3631 bool removeObject, bool removeTool)
3632 {
3633 return booleanOperator(tag, OCC_Internals::Difference, objectDimTags,
3634 toolDimTags, outDimTags, outDimTagsMap, removeObject,
3635 removeTool);
3636 }
3637
booleanFragments(int tag,const std::vector<std::pair<int,int>> & objectDimTags,const std::vector<std::pair<int,int>> & toolDimTags,std::vector<std::pair<int,int>> & outDimTags,std::vector<std::vector<std::pair<int,int>>> & outDimTagsMap,bool removeObject,bool removeTool)3638 bool OCC_Internals::booleanFragments(
3639 int tag, const std::vector<std::pair<int, int> > &objectDimTags,
3640 const std::vector<std::pair<int, int> > &toolDimTags,
3641 std::vector<std::pair<int, int> > &outDimTags,
3642 std::vector<std::vector<std::pair<int, int> > > &outDimTagsMap,
3643 bool removeObject, bool removeTool)
3644 {
3645 return booleanOperator(tag, OCC_Internals::Fragments, objectDimTags,
3646 toolDimTags, outDimTags, outDimTagsMap, removeObject,
3647 removeTool);
3648 }
3649
_getMaxDim()3650 int OCC_Internals::_getMaxDim()
3651 {
3652 if(_tagSolid.Extent()) return 3;
3653 if(_tagFace.Extent()) return 2;
3654 if(_tagEdge.Extent()) return 1;
3655 return 0;
3656 }
3657
_getAllDimTags(std::vector<std::pair<int,int>> & dimTags,int dim)3658 void OCC_Internals::_getAllDimTags(std::vector<std::pair<int, int> > &dimTags,
3659 int dim)
3660 {
3661 for(int d = -2; d < 4; d++) {
3662 if(dim != 99 && dim != d) continue;
3663 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp;
3664 switch(d) {
3665 case 0: exp.Initialize(_tagVertex); break;
3666 case 1: exp.Initialize(_tagEdge); break;
3667 case 2: exp.Initialize(_tagFace); break;
3668 case 3: exp.Initialize(_tagSolid); break;
3669 case -1: exp.Initialize(_tagWire); break;
3670 case -2: exp.Initialize(_tagShell); break;
3671 }
3672 for(; exp.More(); exp.Next())
3673 dimTags.push_back(std::pair<int, int>(d, exp.Key()));
3674 }
3675 }
3676
removeAllDuplicates()3677 void OCC_Internals::removeAllDuplicates()
3678 {
3679 std::vector<std::pair<int, int> > objectDimTags, toolDimTags, outDimTags;
3680 std::vector<std::vector<std::pair<int, int> > > outDimTagsMap;
3681 _getAllDimTags(objectDimTags, _getMaxDim());
3682 booleanFragments(-1, objectDimTags, toolDimTags, outDimTags, outDimTagsMap,
3683 true, true);
3684 }
3685
mergeVertices(const std::vector<int> & tags)3686 bool OCC_Internals::mergeVertices(const std::vector<int> &tags)
3687 {
3688 std::vector<std::pair<int, int> > objectDimTags, toolDimTags, outDimTags;
3689 std::vector<std::vector<std::pair<int, int> > > outDimTagsMap;
3690 for(std::size_t i = 0; i < tags.size(); i++)
3691 objectDimTags.push_back(std::pair<int, int>(0, tags[i]));
3692 return booleanFragments(-1, objectDimTags, toolDimTags, outDimTags,
3693 outDimTagsMap, true, true);
3694 }
3695
_addSimpleShapes(const TopoDS_Shape & shape,std::vector<TopoDS_Shape> & simple)3696 void _addSimpleShapes(const TopoDS_Shape &shape,
3697 std::vector<TopoDS_Shape> &simple)
3698 {
3699 if(shape.ShapeType() != TopAbs_COMPOUND &&
3700 shape.ShapeType() != TopAbs_COMPSOLID) {
3701 simple.push_back(shape);
3702 return;
3703 }
3704
3705 TopTools_MapOfShape mapShape;
3706 TopoDS_Iterator It(shape, Standard_True, Standard_True);
3707
3708 for(; It.More(); It.Next()) {
3709 const TopoDS_Shape &s = It.Value();
3710 if(mapShape.Add(s)) {
3711 if(s.ShapeType() == TopAbs_COMPOUND ||
3712 s.ShapeType() == TopAbs_COMPSOLID) {
3713 _addSimpleShapes(s, simple);
3714 }
3715 else {
3716 simple.push_back(s);
3717 }
3718 }
3719 }
3720 }
3721
_transform(const std::vector<std::pair<int,int>> & inDimTags,BRepBuilderAPI_Transform * tfo,BRepBuilderAPI_GTransform * gtfo)3722 bool OCC_Internals::_transform(
3723 const std::vector<std::pair<int, int> > &inDimTags,
3724 BRepBuilderAPI_Transform *tfo, BRepBuilderAPI_GTransform *gtfo)
3725 {
3726 // build a single compound shape, so that we won't duplicate internal
3727 // boundaries
3728 BRep_Builder b;
3729 TopoDS_Compound c;
3730 b.MakeCompound(c);
3731 for(std::size_t i = 0; i < inDimTags.size(); i++) {
3732 int dim = inDimTags[i].first;
3733 int tag = inDimTags[i].second;
3734 if(!_isBound(dim, tag)) {
3735 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
3736 tag);
3737 return false;
3738 }
3739 TopoDS_Shape shape = _find(dim, tag);
3740 b.Add(c, shape);
3741 }
3742
3743 std::vector<TopoDS_Shape> inShapes;
3744 _addSimpleShapes(c, inShapes);
3745
3746 TopoDS_Shape result;
3747 if(tfo) {
3748 tfo->Perform(c, Standard_False);
3749 if(!tfo->IsDone()) {
3750 Msg::Error("Could not apply transformation");
3751 return false;
3752 }
3753 result = tfo->Shape();
3754 }
3755 else if(gtfo) {
3756 gtfo->Perform(c, Standard_False);
3757 if(!gtfo->IsDone()) {
3758 Msg::Error("Could not apply transformation");
3759 return false;
3760 }
3761 result = gtfo->Shape();
3762 }
3763
3764 // copy vertex-based meshing attributes
3765 TopExp_Explorer exp0;
3766 for(exp0.Init(c, TopAbs_VERTEX); exp0.More(); exp0.Next()) {
3767 TopoDS_Vertex vertex = TopoDS::Vertex(exp0.Current());
3768 TopoDS_Shape transformed;
3769 if(tfo)
3770 transformed = tfo->ModifiedShape(vertex);
3771 else if(gtfo)
3772 transformed = gtfo->ModifiedShape(vertex);
3773 if(!transformed.IsNull()) {
3774 double lc = _attributes->getMeshSize(0, vertex);
3775 if(lc > 0 && lc < MAX_LC)
3776 _attributes->insert(new OCCAttributes(0, transformed, lc));
3777 }
3778 }
3779
3780 // try to re-bind trasnformed shapes to same tags as original shapes
3781 std::vector<TopoDS_Shape> outShapes;
3782 _addSimpleShapes(result, outShapes);
3783
3784 if(inShapes.size() != inDimTags.size() ||
3785 inShapes.size() != outShapes.size()) {
3786 Msg::Error("OpenCASCADE transform changed the number of shapes");
3787 return false;
3788 }
3789 for(std::size_t i = 0; i < inDimTags.size(); i++) {
3790 int dim = inDimTags[i].first;
3791 int tag = inDimTags[i].second;
3792 #if 0
3793 // safe, but slow: _unbind() has linear complexity with respect to the number
3794 // of entities in the model (due to the dependency checking of upward
3795 // adjencencies and the maximum tag update). Using this in a for loop to
3796 // translate copies of entities leads to quadratic complexity.
3797 _unbind(inShapes[i], dim, tag, true);
3798 #else
3799 // bypass it by unbinding the shape and all its subshapes without checking
3800 // dependencies: this is a bit dangerous, as one could translate e.g. the
3801 // face of a cube (this is not allowed!) - which will unbind the face of the
3802 // cube. But the original face will actually be re-bound (with a warning) at
3803 // the next syncronization point, so it's not too bad...
3804 _unbindWithoutChecks(inShapes[i]);
3805 for(int d = -2; d <= 3; d++) _recomputeMaxTag(d);
3806
3807 // FIXME: it would be even better to code a rebind() function to reuse the
3808 // tags not only of the shape, but of all the sub-shapes as well
3809 #endif
3810 _bind(outShapes[i], dim, tag, true);
3811 }
3812
3813 return true;
3814 }
3815
translate(const std::vector<std::pair<int,int>> & inDimTags,double dx,double dy,double dz)3816 bool OCC_Internals::translate(
3817 const std::vector<std::pair<int, int> > &inDimTags, double dx, double dy,
3818 double dz)
3819 {
3820 try {
3821 gp_Trsf t;
3822 t.SetTranslation(gp_Pnt(0, 0, 0), gp_Pnt(dx, dy, dz));
3823 BRepBuilderAPI_Transform tfo(t);
3824 return _transform(inDimTags, &tfo, nullptr);
3825 } catch(Standard_Failure &err) {
3826 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
3827 return false;
3828 }
3829 }
3830
rotate(const std::vector<std::pair<int,int>> & inDimTags,double x,double y,double z,double ax,double ay,double az,double angle)3831 bool OCC_Internals::rotate(const std::vector<std::pair<int, int> > &inDimTags,
3832 double x, double y, double z, double ax, double ay,
3833 double az, double angle)
3834 {
3835 try {
3836 gp_Trsf t;
3837 gp_Ax1 axisOfRevolution(gp_Pnt(x, y, z), gp_Dir(ax, ay, az));
3838 t.SetRotation(axisOfRevolution, angle);
3839 BRepBuilderAPI_Transform tfo(t);
3840 return _transform(inDimTags, &tfo, nullptr);
3841 } catch(Standard_Failure &err) {
3842 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
3843 return false;
3844 }
3845 }
3846
dilate(const std::vector<std::pair<int,int>> & inDimTags,double x,double y,double z,double a,double b,double c)3847 bool OCC_Internals::dilate(const std::vector<std::pair<int, int> > &inDimTags,
3848 double x, double y, double z, double a, double b,
3849 double c)
3850 {
3851 try {
3852 gp_GTrsf gt;
3853 gt.SetVectorialPart(gp_Mat(a, 0, 0, 0, b, 0, 0, 0, c));
3854 gt.SetTranslationPart(gp_XYZ(x * (1 - a), y * (1 - b), z * (1 - c)));
3855 BRepBuilderAPI_GTransform gtfo(gt);
3856 return _transform(inDimTags, nullptr, >fo);
3857 } catch(Standard_Failure &err) {
3858 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
3859 return false;
3860 }
3861 }
3862
symmetry(const std::vector<std::pair<int,int>> & inDimTags,double a,double b,double c,double d)3863 bool OCC_Internals::symmetry(const std::vector<std::pair<int, int> > &inDimTags,
3864 double a, double b, double c, double d)
3865 {
3866 try {
3867 gp_GTrsf gt;
3868 double p = (a * a + b * b + c * c);
3869 if(!p) p = 1e-12;
3870 double f = -2.0 / p;
3871 gt.SetVectorialPart(gp_Mat(1 + a * a * f, a * b * f, a * c * f, a * b * f,
3872 1. + b * b * f, b * c * f, a * c * f, b * c * f,
3873 1. + c * c * f));
3874 gt.SetTranslationPart(gp_XYZ(a * d * f, b * d * f, c * d * f));
3875 BRepBuilderAPI_GTransform gtfo(gt);
3876 return _transform(inDimTags, nullptr, >fo);
3877 } catch(Standard_Failure &err) {
3878 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
3879 return false;
3880 }
3881 }
3882
affine(const std::vector<std::pair<int,int>> & inDimTags,const std::vector<double> & mat)3883 bool OCC_Internals::affine(const std::vector<std::pair<int, int> > &inDimTags,
3884 const std::vector<double> &mat)
3885 {
3886 try {
3887 std::vector<double> a(mat);
3888 if(a.size() < 12) {
3889 Msg::Warning("%d < 12 entries in affine transform matrix", (int)a.size());
3890 a.resize(12, 0.);
3891 }
3892 gp_GTrsf gt;
3893 gt.SetVectorialPart(
3894 gp_Mat(a[0], a[1], a[2], a[4], a[5], a[6], a[8], a[9], a[10]));
3895 gt.SetTranslationPart(gp_XYZ(a[3], a[7], a[11]));
3896 BRepBuilderAPI_GTransform gtfo(gt);
3897 return _transform(inDimTags, nullptr, >fo);
3898 } catch(Standard_Failure &err) {
3899 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
3900 return false;
3901 }
3902 }
3903
copy(const std::vector<std::pair<int,int>> & inDimTags,std::vector<std::pair<int,int>> & outDimTags)3904 bool OCC_Internals::copy(const std::vector<std::pair<int, int> > &inDimTags,
3905 std::vector<std::pair<int, int> > &outDimTags)
3906 {
3907 bool ret = true;
3908 for(std::size_t i = 0; i < inDimTags.size(); i++) {
3909 int dim = inDimTags[i].first;
3910 int tag = inDimTags[i].second;
3911 if(!_isBound(dim, tag)) {
3912 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
3913 tag);
3914 ret = false;
3915 continue;
3916 }
3917 TopoDS_Shape result = BRepBuilderAPI_Copy(_find(dim, tag)).Shape();
3918 int newtag = getMaxTag(dim) + 1;
3919 _bind(result, dim, newtag, true);
3920 outDimTags.push_back(std::make_pair(dim, newtag));
3921 }
3922 return ret;
3923 }
3924
remove(int dim,int tag,bool recursive)3925 bool OCC_Internals::remove(int dim, int tag, bool recursive)
3926 {
3927 if(!_isBound(dim, tag)) {
3928 Msg::Warning("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
3929 tag);
3930 return false;
3931 }
3932 _unbind(_find(dim, tag), dim, tag, recursive);
3933 return true;
3934 }
3935
remove(const std::vector<std::pair<int,int>> & dimTags,bool recursive)3936 bool OCC_Internals::remove(const std::vector<std::pair<int, int> > &dimTags,
3937 bool recursive)
3938 {
3939 bool ret = true;
3940 for(std::size_t i = 0; i < dimTags.size(); i++) {
3941 if(!remove(dimTags[i].first, dimTags[i].second, recursive)) ret = false;
3942 }
3943 return ret;
3944 }
3945
setTargetUnit(const std::string & unit)3946 static void setTargetUnit(const std::string &unit)
3947 {
3948 if(unit.empty()) return; // use unit specified in the file
3949 if(!Interface_Static::SetCVal("xstep.cascade.unit", unit.c_str()))
3950 Msg::Error("Could not set OpenCASCADE target unit '%s'", unit.c_str());
3951 }
3952
3953 #if defined(HAVE_OCC_CAF)
3954
setShapeAttributes(OCCAttributesRTree * attributes,const Handle_XCAFDoc_ShapeTool & shapeTool,const Handle_XCAFDoc_ColorTool & colorTool,const Handle_XCAFDoc_MaterialTool & materialTool,const TDF_Label & label,const TopLoc_Location & loc,const std::string & pathName,bool isRef)3955 static void setShapeAttributes(OCCAttributesRTree *attributes,
3956 const Handle_XCAFDoc_ShapeTool &shapeTool,
3957 const Handle_XCAFDoc_ColorTool &colorTool,
3958 const Handle_XCAFDoc_MaterialTool &materialTool,
3959 const TDF_Label &label,
3960 const TopLoc_Location &loc,
3961 const std::string &pathName, bool isRef)
3962 {
3963 std::string phys = pathName;
3964 Handle(TDataStd_Name) n;
3965 if(label.FindAttribute(TDataStd_Name::GetID(), n)) {
3966 TCollection_ExtendedString name = n->Get();
3967 if(!phys.empty()) phys += "/";
3968 phys += TCollection_AsciiString(name).ToCString();
3969 }
3970
3971 TopLoc_Location partLoc = loc;
3972 Handle(XCAFDoc_Location) l;
3973 if(label.FindAttribute(XCAFDoc_Location::GetID(), l)) {
3974 if(isRef)
3975 partLoc = partLoc * l->Get();
3976 else
3977 partLoc = l->Get();
3978 }
3979
3980 TDF_Label ref;
3981 if(shapeTool->IsReference(label) && shapeTool->GetReferredShape(label, ref)) {
3982 setShapeAttributes(attributes, shapeTool, colorTool, materialTool, ref,
3983 partLoc, phys, true);
3984 }
3985
3986 if(shapeTool->IsSimpleShape(label) && (isRef || shapeTool->IsFree(label))) {
3987 TopoDS_Shape shape = shapeTool->GetShape(label);
3988 shape.Location(isRef ? loc : partLoc);
3989 int dim =
3990 (shape.ShapeType() == TopAbs_VERTEX) ? 0 :
3991 (shape.ShapeType() == TopAbs_EDGE || shape.ShapeType() == TopAbs_WIRE) ?
3992 1 :
3993 (shape.ShapeType() == TopAbs_FACE || shape.ShapeType() == TopAbs_SHELL) ?
3994 2 :
3995 3;
3996
3997 Handle(TCollection_HAsciiString) matName;
3998 Handle(TCollection_HAsciiString) matDescription;
3999 Standard_Real matDensity;
4000 Handle(TCollection_HAsciiString) matDensName;
4001 Handle(TCollection_HAsciiString) matDensValType;
4002 if(materialTool->GetMaterial(label, matName, matDescription, matDensity,
4003 matDensName, matDensValType)) {
4004 if(!phys.empty()) phys += " & ";
4005 phys += matName->ToCString();
4006 Msg::Info(" - Label & material '%s' (%dD)", phys.c_str());
4007 }
4008 else if(phys.size()) {
4009 Msg::Info(" - Label '%s' (%dD)", phys.c_str(), dim);
4010 }
4011 if(phys.size()) { attributes->insert(new OCCAttributes(dim, shape, phys)); }
4012
4013 Quantity_Color col;
4014 if(colorTool->GetColor(label, XCAFDoc_ColorGen, col)) {
4015 double r = col.Red(), g = col.Green(), b = col.Blue();
4016 Msg::Info(" - Color (%g, %g, %g) (%dD)", r, g, b, dim);
4017 attributes->insert(new OCCAttributes(dim, shape, r, g, b, 1.));
4018 }
4019 else if(colorTool->GetColor(label, XCAFDoc_ColorSurf, col)) {
4020 double r = col.Red(), g = col.Green(), b = col.Blue();
4021 Msg::Info(" - Color (%g, %g, %g) (%dD & Surfaces)", r, g, b, dim);
4022 attributes->insert(new OCCAttributes(dim, shape, r, g, b, 1., 1));
4023 }
4024 else if(colorTool->GetColor(label, XCAFDoc_ColorCurv, col)) {
4025 double r = col.Red(), g = col.Green(), b = col.Blue();
4026 Msg::Info(" - Color (%g, %g, %g) (%dD & Curves)", r, g, b, dim);
4027 attributes->insert(new OCCAttributes(dim, shape, r, g, b, 1., 2));
4028 }
4029 // check explicit coloring of boundary entities
4030 if(dim == 3) {
4031 TopExp_Explorer xp2(shape, TopAbs_FACE);
4032 while(xp2.More()) {
4033 if(colorTool->GetColor(xp2.Current(), XCAFDoc_ColorGen, col) ||
4034 colorTool->GetColor(xp2.Current(), XCAFDoc_ColorSurf, col) ||
4035 colorTool->GetColor(xp2.Current(), XCAFDoc_ColorCurv, col)) {
4036 double r = col.Red(), g = col.Green(), b = col.Blue();
4037 Msg::Info(" - Color (%g, %g, %g) (Surface)", r, g, b);
4038 TopoDS_Face face = TopoDS::Face(xp2.Current());
4039 attributes->insert(new OCCAttributes(2, face, r, g, b, 1.));
4040 }
4041 xp2.Next();
4042 }
4043 }
4044 if(dim == 2) {
4045 TopExp_Explorer xp1(shape, TopAbs_EDGE);
4046 while(xp1.More()) {
4047 if(colorTool->GetColor(xp1.Current(), XCAFDoc_ColorGen, col) ||
4048 colorTool->GetColor(xp1.Current(), XCAFDoc_ColorSurf, col) ||
4049 colorTool->GetColor(xp1.Current(), XCAFDoc_ColorCurv, col)) {
4050 double r = col.Red(), g = col.Green(), b = col.Blue();
4051 Msg::Info(" - Color (%g, %g, %g) (Curve)", r, g, b);
4052 attributes->insert(
4053 new OCCAttributes(1, TopoDS::Face(xp1.Current()), r, g, b, 1.));
4054 }
4055 xp1.Next();
4056 }
4057 }
4058 }
4059 else {
4060 for(TDF_ChildIterator it(label); it.More(); it.Next()) {
4061 setShapeAttributes(attributes, shapeTool, colorTool, materialTool,
4062 it.Value(), partLoc, phys, isRef);
4063 }
4064 }
4065 }
4066
4067 template <class T>
readAttributes(OCCAttributesRTree * attributes,T & reader,const std::string & format)4068 void readAttributes(OCCAttributesRTree *attributes, T &reader,
4069 const std::string &format)
4070 {
4071 // dummy XCAF Application to handle the STEP XCAF Document
4072 static Handle_XCAFApp_Application dummy_app =
4073 XCAFApp_Application::GetApplication();
4074 // XCAF Document to contain the STEP/IGES file itself
4075 Handle_TDocStd_Document doc;
4076 // check if a file is already open under this handle, if so, close it to
4077 // prevent segfaults when trying to create a new document
4078 if(dummy_app->NbDocuments() > 0) {
4079 dummy_app->GetDocument(1, doc);
4080 dummy_app->Close(doc);
4081 }
4082 dummy_app->NewDocument(format.c_str(), doc);
4083 // transfer STEP/IGES into the document, and get the main label
4084 reader.Transfer(doc);
4085 TDF_Label mainLabel = doc->Main();
4086 Handle_XCAFDoc_ShapeTool shapeTool =
4087 XCAFDoc_DocumentTool::ShapeTool(mainLabel);
4088 Handle_XCAFDoc_ColorTool colorTool =
4089 XCAFDoc_DocumentTool::ColorTool(mainLabel);
4090 Handle_XCAFDoc_MaterialTool materialTool =
4091 XCAFDoc_DocumentTool::MaterialTool(mainLabel);
4092 // traverse the labels recursively to set attributes on shapes
4093 setShapeAttributes(attributes, shapeTool, colorTool, materialTool, mainLabel,
4094 TopLoc_Location(), "", false);
4095 }
4096
4097 #endif
4098
importShapes(const std::string & fileName,bool highestDimOnly,std::vector<std::pair<int,int>> & outDimTags,const std::string & format)4099 bool OCC_Internals::importShapes(const std::string &fileName,
4100 bool highestDimOnly,
4101 std::vector<std::pair<int, int> > &outDimTags,
4102 const std::string &format)
4103 {
4104 std::vector<std::string> split = SplitFileName(fileName);
4105
4106 TCollection_AsciiString occfile(fileName.c_str());
4107
4108 TopoDS_Shape result;
4109 try {
4110 if(format == "brep" || split[2] == ".brep" || split[2] == ".BREP") {
4111 BRep_Builder aBuilder;
4112 BRepTools::Read(result, occfile.ToCString(), aBuilder);
4113 }
4114 else if(format == "step" || split[2] == ".step" || split[2] == ".stp" ||
4115 split[2] == ".STEP" || split[2] == ".STP") {
4116 STEPControl_Reader reader;
4117 setTargetUnit(CTX::instance()->geom.occTargetUnit);
4118 #if defined(HAVE_OCC_CAF)
4119 //Interface_Static::SetIVal("read.stepcaf.subshapes.name", 1);
4120 STEPCAFControl_Reader cafreader;
4121 if(cafreader.ReadFile(occfile.ToCString()) != IFSelect_RetDone) {
4122 Msg::Error("Could not read file '%s'", fileName.c_str());
4123 return false;
4124 }
4125 if(CTX::instance()->geom.occImportLabels)
4126 readAttributes(_attributes, cafreader, "STEP-XCAF");
4127 reader = cafreader.ChangeReader();
4128 #else
4129 if(reader.ReadFile(occfile.ToCString()) != IFSelect_RetDone) {
4130 Msg::Error("Could not read file '%s'", fileName.c_str());
4131 return false;
4132 }
4133 #endif
4134 reader.NbRootsForTransfer();
4135 reader.TransferRoots();
4136 result = reader.OneShape();
4137 }
4138 else if(format == "iges" || split[2] == ".iges" || split[2] == ".igs" ||
4139 split[2] == ".IGES" || split[2] == ".IGS") {
4140 setTargetUnit(CTX::instance()->geom.occTargetUnit);
4141 #if defined(HAVE_OCC_CAF)
4142 IGESCAFControl_Reader reader;
4143 if(reader.ReadFile(occfile.ToCString()) != IFSelect_RetDone) {
4144 Msg::Error("Could not read file '%s'", fileName.c_str());
4145 return false;
4146 }
4147 if(CTX::instance()->geom.occImportLabels)
4148 readAttributes(_attributes, reader, "IGES-XCAF");
4149 #else
4150 IGESControl_Reader reader;
4151 if(reader.ReadFile(occfile.ToCString()) != IFSelect_RetDone) {
4152 Msg::Error("Could not read file '%s'", fileName.c_str());
4153 return false;
4154 }
4155 #endif
4156 reader.NbRootsForTransfer();
4157 reader.TransferRoots();
4158 result = reader.OneShape();
4159 }
4160 else {
4161 Msg::Error("Unknown file type '%s'", fileName.c_str());
4162 return false;
4163 }
4164 } catch(Standard_Failure &err) {
4165 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
4166 return false;
4167 }
4168
4169 BRepTools::Clean(result);
4170
4171 _healShape(
4172 result, CTX::instance()->geom.tolerance,
4173 CTX::instance()->geom.occFixDegenerated,
4174 CTX::instance()->geom.occFixSmallEdges,
4175 CTX::instance()->geom.occFixSmallFaces, CTX::instance()->geom.occSewFaces,
4176 CTX::instance()->geom.occMakeSolids, CTX::instance()->geom.occScaling);
4177
4178 _multiBind(result, -1, outDimTags, highestDimOnly, true);
4179 return true;
4180 }
4181
importShapes(const TopoDS_Shape * shape,bool highestDimOnly,std::vector<std::pair<int,int>> & outDimTags)4182 bool OCC_Internals::importShapes(const TopoDS_Shape *shape, bool highestDimOnly,
4183 std::vector<std::pair<int, int> > &outDimTags)
4184 {
4185 if(!shape) return false;
4186 _multiBind(*shape, -1, outDimTags, highestDimOnly, true);
4187 return true;
4188 }
4189
exportShapes(GModel * model,const std::string & fileName,const std::string & format)4190 bool OCC_Internals::exportShapes(GModel *model, const std::string &fileName,
4191 const std::string &format)
4192 {
4193 // put all top-level OCC shapes from a GModel in a single compound (we use the
4194 // topology to only consider top-level shapes; otherwise we get duplicates in
4195 // the STEP export)
4196 BRep_Builder b;
4197 TopoDS_Compound c;
4198 b.MakeCompound(c);
4199 for(auto it = model->firstRegion(); it != model->lastRegion(); it++) {
4200 GRegion *gr = *it;
4201 if(gr->getNativeType() == GEntity::OpenCascadeModel)
4202 b.Add(c, *(TopoDS_Solid *)gr->getNativePtr());
4203 }
4204 for(auto it = model->firstFace(); it != model->lastFace(); it++) {
4205 GFace *gf = *it;
4206 if(!gf->numRegions() && gf->getNativeType() == GEntity::OpenCascadeModel)
4207 b.Add(c, *(TopoDS_Face *)gf->getNativePtr());
4208 }
4209 for(auto it = model->firstEdge(); it != model->lastEdge(); it++) {
4210 GEdge *ge = *it;
4211 if(!ge->numFaces() && ge->getNativeType() == GEntity::OpenCascadeModel)
4212 b.Add(c, *(TopoDS_Edge *)ge->getNativePtr());
4213 }
4214 for(auto it = model->firstVertex(); it != model->lastVertex(); it++) {
4215 GVertex *gv = *it;
4216 if(!gv->numEdges() && gv->getNativeType() == GEntity::OpenCascadeModel)
4217 b.Add(c, *(TopoDS_Vertex *)gv->getNativePtr());
4218 }
4219
4220 std::vector<std::string> split = SplitFileName(fileName);
4221
4222 TCollection_AsciiString occfile(fileName.c_str());
4223
4224 try {
4225 if(format == "brep" || split[2] == ".brep" || split[2] == ".BREP") {
4226 BRepTools::Write(c, occfile.ToCString());
4227 }
4228 else if(format == "step" || split[2] == ".step" || split[2] == ".stp" ||
4229 split[2] == ".STEP" || split[2] == ".STP") {
4230 STEPControl_Writer writer;
4231 setTargetUnit(CTX::instance()->geom.occTargetUnit);
4232 if(writer.Transfer(c, STEPControl_AsIs) == IFSelect_RetDone) {
4233 if(writer.Write(occfile.ToCString()) != IFSelect_RetDone) {
4234 Msg::Error("Could not create file '%s'", fileName.c_str());
4235 return false;
4236 }
4237 }
4238 else {
4239 Msg::Error("Could not create STEP data");
4240 return false;
4241 }
4242 }
4243 } catch(Standard_Failure &err) {
4244 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
4245 return false;
4246 }
4247 return true;
4248 }
4249
setMeshSize(int dim,int tag,double size)4250 void OCC_Internals::setMeshSize(int dim, int tag, double size)
4251 {
4252 if(dim != 0) return;
4253 if(_tagVertex.IsBound(tag)) {
4254 OCCAttributes *a = new OCCAttributes(0, _tagVertex.Find(tag), size);
4255 // first remove any other constraint
4256 _attributes->remove(a);
4257 _attributes->insert(a);
4258 }
4259 }
4260
getEntities(std::vector<std::pair<int,int>> & dimTags,int dim)4261 bool OCC_Internals::getEntities(std::vector<std::pair<int, int> > &dimTags,
4262 int dim)
4263 {
4264 for(int d = 0; d < 4; d++) {
4265 if(dim != -1 && dim != d) continue;
4266 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp;
4267 switch(d) {
4268 case 0: exp.Initialize(_tagVertex); break;
4269 case 1: exp.Initialize(_tagEdge); break;
4270 case 2: exp.Initialize(_tagFace); break;
4271 case 3: exp.Initialize(_tagSolid); break;
4272 }
4273 for(; exp.More(); exp.Next())
4274 dimTags.push_back(std::make_pair(d, exp.Key()));
4275 }
4276 return true;
4277 }
4278
getVertex(int tag,double & x,double & y,double & z)4279 bool OCC_Internals::getVertex(int tag, double &x, double &y, double &z)
4280 {
4281 if(_tagVertex.IsBound(tag)) {
4282 gp_Pnt pnt = BRep_Tool::Pnt(TopoDS::Vertex(_tagVertex.Find(tag)));
4283 x = pnt.X();
4284 y = pnt.Y();
4285 z = pnt.Z();
4286 return true;
4287 }
4288 return false;
4289 }
4290
_getBoundingBox(const TopoDS_Shape & shape,double & xmin,double & ymin,double & zmin,double & xmax,double & ymax,double & zmax)4291 bool OCC_Internals::_getBoundingBox(const TopoDS_Shape &shape, double &xmin,
4292 double &ymin, double &zmin, double &xmax,
4293 double &ymax, double &zmax)
4294 {
4295 if(CTX::instance()->geom.occBoundsUseSTL) {
4296 std::vector<SPoint3> vertices;
4297 std::vector<SVector3> normals;
4298 std::vector<int> triangles;
4299 _makeSTL(shape, vertices, normals, triangles);
4300 }
4301 Bnd_Box b;
4302 try {
4303 BRepBndLib::Add(shape, b);
4304 } catch(Standard_Failure &err) {
4305 Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
4306 return false;
4307 }
4308 b.Get(xmin, ymin, zmin, xmax, ymax, zmax);
4309 if(CTX::instance()->geom.occBoundsUseSTL)
4310 fixSTLBounds(xmin, ymin, zmin, xmax, ymax, zmax);
4311 return true;
4312 }
4313
getBoundingBox(int dim,int tag,double & xmin,double & ymin,double & zmin,double & xmax,double & ymax,double & zmax)4314 bool OCC_Internals::getBoundingBox(int dim, int tag, double &xmin, double &ymin,
4315 double &zmin, double &xmax, double &ymax,
4316 double &zmax)
4317 {
4318 if(!_isBound(dim, tag)) {
4319 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
4320 tag);
4321 return false;
4322 }
4323 TopoDS_Shape shape = _find(dim, tag);
4324 return _getBoundingBox(shape, xmin, ymin, zmin, xmax, ymax, zmax);
4325 }
4326
getEntitiesInBoundingBox(double xmin,double ymin,double zmin,double xmax,double ymax,double zmax,std::vector<std::pair<int,int>> & dimTags,int dim)4327 bool OCC_Internals::getEntitiesInBoundingBox(
4328 double xmin, double ymin, double zmin, double xmax, double ymax, double zmax,
4329 std::vector<std::pair<int, int> > &dimTags, int dim)
4330 {
4331 // if we use this often, create an rtree to avoid the linear search
4332 for(int d = 0; d < 4; d++) {
4333 if(dim != -1 && dim != d) continue;
4334 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp;
4335 switch(d) {
4336 case 0: exp.Initialize(_tagVertex); break;
4337 case 1: exp.Initialize(_tagEdge); break;
4338 case 2: exp.Initialize(_tagFace); break;
4339 case 3: exp.Initialize(_tagSolid); break;
4340 }
4341 for(; exp.More(); exp.Next()) {
4342 double xmin2 = 0, ymin2 = 0, zmin2 = 0, xmax2 = 0, ymax2 = 0, zmax2 = 0;
4343 _getBoundingBox(exp.Value(), xmin2, ymin2, zmin2, xmax2, ymax2, zmax2);
4344 if(xmin2 >= xmin && xmax2 <= xmax && ymin2 >= ymin && ymax2 <= ymax &&
4345 zmin2 >= zmin && zmax2 <= zmax)
4346 dimTags.push_back(std::make_pair(d, exp.Key()));
4347 }
4348 }
4349 return true;
4350 }
4351
getMass(int dim,int tag,double & mass)4352 bool OCC_Internals::getMass(int dim, int tag, double &mass)
4353 {
4354 if(!_isBound(dim, tag)) {
4355 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
4356 tag);
4357 return false;
4358 }
4359 TopoDS_Shape shape = _find(dim, tag);
4360 GProp_GProps System;
4361 switch(dim) {
4362 case 1: BRepGProp::LinearProperties(shape, System); break;
4363 case 2: BRepGProp::SurfaceProperties(shape, System); break;
4364 case 3: BRepGProp::VolumeProperties(shape, System); break;
4365 }
4366 mass = System.Mass();
4367 return true;
4368 }
4369
getCenterOfMass(int dim,int tag,double & x,double & y,double & z)4370 bool OCC_Internals::getCenterOfMass(int dim, int tag, double &x, double &y,
4371 double &z)
4372 {
4373 if(!_isBound(dim, tag)) {
4374 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
4375 tag);
4376 return false;
4377 }
4378 TopoDS_Shape shape = _find(dim, tag);
4379 GProp_GProps System;
4380 switch(dim) {
4381 case 1: BRepGProp::LinearProperties(shape, System); break;
4382 case 2: BRepGProp::SurfaceProperties(shape, System); break;
4383 case 3: BRepGProp::VolumeProperties(shape, System); break;
4384 }
4385 gp_Pnt c = System.CentreOfMass();
4386 x = c.X();
4387 y = c.Y();
4388 z = c.Z();
4389 return true;
4390 }
4391
getMatrixOfInertia(int dim,int tag,std::vector<double> & mat)4392 bool OCC_Internals::getMatrixOfInertia(int dim, int tag,
4393 std::vector<double> &mat)
4394 {
4395 if(!_isBound(dim, tag)) {
4396 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d", dim,
4397 tag);
4398 return false;
4399 }
4400 TopoDS_Shape shape = _find(dim, tag);
4401 GProp_GProps System;
4402 switch(dim) {
4403 case 1: BRepGProp::LinearProperties(shape, System); break;
4404 case 2: BRepGProp::SurfaceProperties(shape, System); break;
4405 case 3: BRepGProp::VolumeProperties(shape, System); break;
4406 }
4407 gp_Mat m = System.MatrixOfInertia();
4408 mat.clear();
4409 for(int i = 1; i <= 3; i++)
4410 for(int j = 1; j <= 3; j++) mat.push_back(m.Value(i, j));
4411 return true;
4412 }
4413
getDistance(int dim1,int tag1,int dim2,int tag2,double & x1,double & y1,double & z1,double & x2,double & y2,double & z2)4414 double OCC_Internals::getDistance(int dim1, int tag1,
4415 int dim2, int tag2,
4416 double &x1, double &y1, double &z1,
4417 double &x2, double &y2, double &z2)
4418 {
4419 if(!_isBound(dim1, tag1)) {
4420 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d",
4421 dim1, tag1);
4422 return false;
4423 }
4424 TopoDS_Shape shape1 = _find(dim1, tag1);
4425
4426 if(!_isBound(dim2, tag2)) {
4427 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d",
4428 dim2, tag2);
4429 return false;
4430 }
4431 TopoDS_Shape shape2 = _find(dim2, tag2);
4432
4433 BRepExtrema_DistShapeShape dist(shape1, shape2);
4434 if(dist.IsDone()) {
4435 double dmin = 1.e200;
4436 gp_Pnt pmin1, pmin2;
4437 for(int i = 1; i <= dist.NbSolution(); i++) {
4438 gp_Pnt p1 = dist.PointOnShape1(i);
4439 gp_Pnt p2 = dist.PointOnShape2(i);
4440 double d = p1.Distance(p2);
4441 if(d < dmin) {
4442 dmin = d;
4443 pmin1 = p1;
4444 pmin2 = p2;
4445 }
4446 }
4447 x1 = pmin1.X();
4448 y1 = pmin1.Y();
4449 z1 = pmin1.Z();
4450 x2 = pmin2.X();
4451 y2 = pmin2.Y();
4452 z2 = pmin2.Z();
4453 return dmin;
4454 }
4455
4456 return -1.;
4457 }
4458
sortByInvDim(std::pair<int,int> const & lhs,std::pair<int,int> const & rhs)4459 bool const sortByInvDim(std::pair<int, int> const &lhs,
4460 std::pair<int, int> const &rhs)
4461 {
4462 return lhs.first > rhs.first;
4463 }
4464
synchronize(GModel * model)4465 void OCC_Internals::synchronize(GModel *model)
4466 {
4467 Msg::Debug("Syncing OCC_Internals with GModel");
4468
4469 // make sure to remove from GModel all entities that have been deleted in
4470 // OCC_Internals since the last synchronization
4471 std::vector<std::pair<int, int> > toRemove;
4472 toRemove.insert(toRemove.end(), _toRemove.begin(), _toRemove.end());
4473 Msg::Debug("Sync is removing %d model entities", toRemove.size());
4474 // make sure to delete highest dimensional entities first (model->remove()
4475 // will not remove entities that are the boundary of others!)
4476 std::sort(toRemove.begin(), toRemove.end(), sortByInvDim);
4477 std::vector<GEntity*> removed;
4478 model->remove(toRemove, removed);
4479 Msg::Debug("Destroying %lu entities in model", removed.size());
4480 for(std::size_t i = 0; i < removed.size(); i++) delete removed[i];
4481 _toRemove.clear();
4482
4483 // iterate over all shapes with tags, and import them into the (sub)shape
4484 // _maps
4485 _somap.Clear();
4486 _shmap.Clear();
4487 _fmap.Clear();
4488 _wmap.Clear();
4489 _emap.Clear();
4490 _vmap.Clear();
4491 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp0(_tagVertex);
4492 for(; exp0.More(); exp0.Next()) _addShapeToMaps(exp0.Value());
4493 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp1(_tagEdge);
4494 for(; exp1.More(); exp1.Next()) _addShapeToMaps(exp1.Value());
4495 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp2(_tagFace);
4496 for(; exp2.More(); exp2.Next()) _addShapeToMaps(exp2.Value());
4497 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp3(_tagSolid);
4498 for(; exp3.More(); exp3.Next()) _addShapeToMaps(exp3.Value());
4499
4500 // import all shapes in _maps into the GModel, preserving all explicit tags
4501 int vTagMax = std::max(model->getMaxElementaryNumber(0), getMaxTag(0));
4502 int eTagMax = std::max(model->getMaxElementaryNumber(1), getMaxTag(1));
4503 int fTagMax = std::max(model->getMaxElementaryNumber(2), getMaxTag(2));
4504 int rTagMax = std::max(model->getMaxElementaryNumber(3), getMaxTag(3));
4505 for(int i = 1; i <= _vmap.Extent(); i++) {
4506 TopoDS_Vertex vertex = TopoDS::Vertex(_vmap(i));
4507 GVertex *occv = getVertexForOCCShape(model, vertex);
4508 if(!occv) {
4509 int tag;
4510 if(_vertexTag.IsBound(vertex))
4511 tag = _vertexTag.Find(vertex);
4512 else {
4513 tag = ++vTagMax;
4514 Msg::Debug("Binding unbound OpenCASCADE point to tag %d", tag);
4515 _bind(vertex, tag);
4516 }
4517 occv = new OCCVertex(model, vertex, tag);
4518 model->add(occv);
4519 }
4520 double lc = _attributes->getMeshSize(0, vertex);
4521 if(lc != MAX_LC) occv->setPrescribedMeshSizeAtVertex(lc);
4522 std::vector<std::string> labels;
4523 _attributes->getLabels(0, vertex, labels);
4524 if(labels.size()) model->setElementaryName(0, occv->tag(), labels[0]);
4525 unsigned int col = 0, boundary = 0;
4526 if(_attributes->getColor(0, vertex, col, boundary)) { occv->setColor(col); }
4527 }
4528 for(int i = 1; i <= _emap.Extent(); i++) {
4529 TopoDS_Edge edge = TopoDS::Edge(_emap(i));
4530 GEdge *occe = getEdgeForOCCShape(model, edge);
4531 if(!occe) {
4532 GVertex *v1 = getVertexForOCCShape(model, TopExp::FirstVertex(edge));
4533 GVertex *v2 = getVertexForOCCShape(model, TopExp::LastVertex(edge));
4534 int tag;
4535 if(_edgeTag.IsBound(edge))
4536 tag = _edgeTag.Find(edge);
4537 else {
4538 tag = ++eTagMax;
4539 Msg::Debug("Binding unbound OpenCASCADE curve to tag %d", tag);
4540 _bind(edge, tag);
4541 }
4542 occe = new OCCEdge(model, edge, tag, v1, v2);
4543 model->add(occe);
4544 }
4545 _copyExtrudedAttributes(edge, occe);
4546 std::vector<std::string> labels;
4547 _attributes->getLabels(1, edge, labels);
4548 if(labels.size()) model->setElementaryName(1, occe->tag(), labels[0]);
4549 unsigned int col = 0, boundary = 0;
4550 if(_attributes->getColor(1, edge, col, boundary)) { occe->setColor(col); }
4551 }
4552 for(int i = 1; i <= _fmap.Extent(); i++) {
4553 TopoDS_Face face = TopoDS::Face(_fmap(i));
4554 GFace *occf = getFaceForOCCShape(model, face);
4555 if(!occf) {
4556 int tag;
4557 if(_faceTag.IsBound(face))
4558 tag = _faceTag.Find(face);
4559 else {
4560 tag = ++fTagMax;
4561 Msg::Debug("Binding unbound OpenCASCADE surface to tag %d", tag);
4562 _bind(face, tag);
4563 }
4564 occf = new OCCFace(model, face, tag);
4565 model->add(occf);
4566 }
4567 _copyExtrudedAttributes(face, occf);
4568 std::vector<std::string> labels;
4569 _attributes->getLabels(2, face, labels);
4570 if(labels.size()) model->setElementaryName(2, occf->tag(), labels[0]);
4571 unsigned int col = 0, boundary = 0;
4572 if(_attributes->getColor(2, face, col, boundary)) {
4573 occf->setColor(col);
4574 if(boundary == 2) {
4575 std::vector<GEdge *> edges = occf->edges();
4576 for(std::size_t j = 0; j < edges.size(); j++) {
4577 // only if not specified explicitly before
4578 if(!edges[j]->useColor()) edges[j]->setColor(col);
4579 }
4580 }
4581 }
4582 }
4583 for(int i = 1; i <= _somap.Extent(); i++) {
4584 TopoDS_Solid region = TopoDS::Solid(_somap(i));
4585 GRegion *occr = getRegionForOCCShape(model, region);
4586 if(!occr) {
4587 int tag;
4588 if(_solidTag.IsBound(region))
4589 tag = _solidTag(region);
4590 else {
4591 tag = ++rTagMax;
4592 Msg::Debug("Binding unbound OpenCASCADE volume to tag %d", tag);
4593 _bind(region, tag);
4594 }
4595 occr = new OCCRegion(model, region, tag);
4596 model->add(occr);
4597 }
4598 _copyExtrudedAttributes(region, occr);
4599 std::vector<std::string> labels;
4600 _attributes->getLabels(3, region, labels);
4601 if(labels.size()) model->setElementaryName(3, occr->tag(), labels[0]);
4602 unsigned int col = 0, boundary = 0;
4603 if(_attributes->getColor(3, region, col, boundary)) {
4604 occr->setColor(col);
4605 if(boundary == 1) {
4606 std::vector<GFace *> faces = occr->faces();
4607 for(std::size_t j = 0; j < faces.size(); j++) {
4608 // only if not specified explicitly before
4609 if(!faces[j]->useColor()) faces[j]->setColor(col);
4610 }
4611 }
4612 else if(boundary == 2) {
4613 std::vector<GEdge *> edges = occr->edges();
4614 for(std::size_t j = 0; j < edges.size(); j++) {
4615 // only if not specified explicitly before
4616 if(!edges[j]->useColor()) edges[j]->setColor(col);
4617 }
4618 }
4619 }
4620 }
4621
4622 // if fuzzy boolean tolerance was used, some vertex positions should be
4623 // recomputed (e.g. end point of curves
4624 if(CTX::instance()->geom.toleranceBoolean) model->snapVertices();
4625
4626 // recompute global boundind box in CTX
4627 SetBoundingBox();
4628
4629 Msg::Debug("GModel imported:");
4630 Msg::Debug("%d points", model->getNumVertices());
4631 Msg::Debug("%d curves", model->getNumEdges());
4632 Msg::Debug("%d surfaces", model->getNumFaces());
4633 Msg::Debug("%d volumes", model->getNumRegions());
4634 _changed = false;
4635 }
4636
getVertexForOCCShape(GModel * model,const TopoDS_Vertex & toFind)4637 GVertex *OCC_Internals::getVertexForOCCShape(GModel *model,
4638 const TopoDS_Vertex &toFind)
4639 {
4640 if(_vertexTag.IsBound(toFind))
4641 return model->getVertexByTag(_vertexTag.Find(toFind));
4642 return nullptr;
4643 }
4644
getEdgeForOCCShape(GModel * model,const TopoDS_Edge & toFind)4645 GEdge *OCC_Internals::getEdgeForOCCShape(GModel *model,
4646 const TopoDS_Edge &toFind)
4647 {
4648 if(_edgeTag.IsBound(toFind))
4649 return model->getEdgeByTag(_edgeTag.Find(toFind));
4650 return nullptr;
4651 }
4652
getFaceForOCCShape(GModel * model,const TopoDS_Face & toFind)4653 GFace *OCC_Internals::getFaceForOCCShape(GModel *model,
4654 const TopoDS_Face &toFind)
4655 {
4656 if(_faceTag.IsBound(toFind))
4657 return model->getFaceByTag(_faceTag.Find(toFind));
4658 return nullptr;
4659 }
4660
getRegionForOCCShape(GModel * model,const TopoDS_Solid & toFind)4661 GRegion *OCC_Internals::getRegionForOCCShape(GModel *model,
4662 const TopoDS_Solid &toFind)
4663 {
4664 if(_solidTag.IsBound(toFind))
4665 return model->getRegionByTag(_solidTag.Find(toFind));
4666 return nullptr;
4667 }
4668
_addShapeToMaps(const TopoDS_Shape & shape)4669 void OCC_Internals::_addShapeToMaps(const TopoDS_Shape &shape)
4670 {
4671 // Solids
4672 TopExp_Explorer exp0, exp1, exp2, exp3, exp4, exp5;
4673 for(exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()) {
4674 TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
4675 if(_somap.FindIndex(solid) < 1) {
4676 _somap.Add(solid);
4677 for(exp1.Init(solid, TopAbs_SHELL); exp1.More(); exp1.Next()) {
4678 TopoDS_Shell shell = TopoDS::Shell(exp1.Current());
4679 if(_shmap.FindIndex(shell) < 1) {
4680 _shmap.Add(shell);
4681
4682 for(exp2.Init(shell, TopAbs_FACE); exp2.More(); exp2.Next()) {
4683 TopoDS_Face face = TopoDS::Face(exp2.Current());
4684 if(_fmap.FindIndex(face) < 1) {
4685 _fmap.Add(face);
4686
4687 for(exp3.Init(face.Oriented(TopAbs_FORWARD), TopAbs_WIRE);
4688 exp3.More(); exp3.Next()) {
4689 // for(exp3.Init(face, TopAbs_WIRE); exp3.More(); exp3.Next()){
4690 TopoDS_Wire wire = TopoDS::Wire(exp3.Current());
4691 if(_wmap.FindIndex(wire) < 1) {
4692 _wmap.Add(wire);
4693
4694 for(exp4.Init(wire, TopAbs_EDGE); exp4.More(); exp4.Next()) {
4695 TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
4696 if(_emap.FindIndex(edge) < 1) {
4697 _emap.Add(edge);
4698
4699 for(exp5.Init(edge, TopAbs_VERTEX); exp5.More();
4700 exp5.Next()) {
4701 TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
4702 if(_vmap.FindIndex(vertex) < 1) _vmap.Add(vertex);
4703 }
4704 }
4705 }
4706 }
4707 }
4708 }
4709 }
4710 }
4711 }
4712 }
4713 }
4714
4715 // Free Shells
4716 for(exp1.Init(shape, TopAbs_SHELL, TopAbs_SOLID); exp1.More(); exp1.Next()) {
4717 const TopoDS_Shape &shell = exp1.Current();
4718 if(_shmap.FindIndex(shell) < 1) {
4719 _shmap.Add(shell);
4720
4721 for(exp2.Init(shell, TopAbs_FACE); exp2.More(); exp2.Next()) {
4722 TopoDS_Face face = TopoDS::Face(exp2.Current());
4723 if(_fmap.FindIndex(face) < 1) {
4724 _fmap.Add(face);
4725
4726 for(exp3.Init(face, TopAbs_WIRE); exp3.More(); exp3.Next()) {
4727 TopoDS_Wire wire = TopoDS::Wire(exp3.Current());
4728 if(_wmap.FindIndex(wire) < 1) {
4729 _wmap.Add(wire);
4730
4731 for(exp4.Init(wire, TopAbs_EDGE); exp4.More(); exp4.Next()) {
4732 TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
4733 if(_emap.FindIndex(edge) < 1) {
4734 _emap.Add(edge);
4735
4736 for(exp5.Init(edge, TopAbs_VERTEX); exp5.More();
4737 exp5.Next()) {
4738 TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
4739 if(_vmap.FindIndex(vertex) < 1) _vmap.Add(vertex);
4740 }
4741 }
4742 }
4743 }
4744 }
4745 }
4746 }
4747 }
4748 }
4749
4750 // Free Faces
4751 for(exp2.Init(shape, TopAbs_FACE, TopAbs_SHELL); exp2.More(); exp2.Next()) {
4752 TopoDS_Face face = TopoDS::Face(exp2.Current());
4753 if(_fmap.FindIndex(face) < 1) {
4754 _fmap.Add(face);
4755
4756 for(exp3.Init(face, TopAbs_WIRE); exp3.More(); exp3.Next()) {
4757 TopoDS_Wire wire = TopoDS::Wire(exp3.Current());
4758 if(_wmap.FindIndex(wire) < 1) {
4759 _wmap.Add(wire);
4760
4761 for(exp4.Init(wire, TopAbs_EDGE); exp4.More(); exp4.Next()) {
4762 TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
4763 if(_emap.FindIndex(edge) < 1) {
4764 _emap.Add(edge);
4765
4766 for(exp5.Init(edge, TopAbs_VERTEX); exp5.More(); exp5.Next()) {
4767 TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
4768 if(_vmap.FindIndex(vertex) < 1) _vmap.Add(vertex);
4769 }
4770 }
4771 }
4772 }
4773 }
4774 }
4775 }
4776
4777 // Free Wires
4778 for(exp3.Init(shape, TopAbs_WIRE, TopAbs_FACE); exp3.More(); exp3.Next()) {
4779 TopoDS_Wire wire = TopoDS::Wire(exp3.Current());
4780 if(_wmap.FindIndex(wire) < 1) {
4781 _wmap.Add(wire);
4782
4783 for(exp4.Init(wire, TopAbs_EDGE); exp4.More(); exp4.Next()) {
4784 TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
4785 if(_emap.FindIndex(edge) < 1) {
4786 _emap.Add(edge);
4787
4788 for(exp5.Init(edge, TopAbs_VERTEX); exp5.More(); exp5.Next()) {
4789 TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
4790 if(_vmap.FindIndex(vertex) < 1) _vmap.Add(vertex);
4791 }
4792 }
4793 }
4794 }
4795 }
4796
4797 // Free Edges
4798 for(exp4.Init(shape, TopAbs_EDGE, TopAbs_WIRE); exp4.More(); exp4.Next()) {
4799 TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
4800 if(_emap.FindIndex(edge) < 1) {
4801 _emap.Add(edge);
4802
4803 for(exp5.Init(edge, TopAbs_VERTEX); exp5.More(); exp5.Next()) {
4804 TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
4805 if(_vmap.FindIndex(vertex) < 1) _vmap.Add(vertex);
4806 }
4807 }
4808 }
4809
4810 // Free Vertices
4811 for(exp5.Init(shape, TopAbs_VERTEX, TopAbs_EDGE); exp5.More(); exp5.Next()) {
4812 TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
4813 if(_vmap.FindIndex(vertex) < 1) { _vmap.Add(vertex); }
4814 }
4815 }
4816
_healShape(TopoDS_Shape & myshape,double tolerance,bool fixDegenerated,bool fixSmallEdges,bool fixSmallFaces,bool sewFaces,bool makeSolids,double scaling)4817 void OCC_Internals::_healShape(TopoDS_Shape &myshape, double tolerance,
4818 bool fixDegenerated, bool fixSmallEdges,
4819 bool fixSmallFaces, bool sewFaces,
4820 bool makeSolids, double scaling)
4821 {
4822 if(scaling != 1.0) {
4823 Msg::Info("Scaling geometry (factor: %g)", scaling);
4824 gp_Trsf t;
4825 t.SetScaleFactor(scaling);
4826 BRepBuilderAPI_Transform trsf(myshape, t);
4827 myshape = trsf.Shape();
4828 }
4829
4830 if(!fixDegenerated && !fixSmallEdges && !fixSmallFaces && !sewFaces &&
4831 !makeSolids)
4832 return;
4833
4834 Msg::Info("Healing shapes (tolerance: %g)", tolerance);
4835 double t1 = Cpu(), w1 = TimeOfDay();
4836
4837 _somap.Clear();
4838 _shmap.Clear();
4839 _fmap.Clear();
4840 _wmap.Clear();
4841 _emap.Clear();
4842 _vmap.Clear();
4843 _addShapeToMaps(myshape);
4844
4845 TopExp_Explorer exp0, exp1;
4846 int nrc = 0, nrcs = 0;
4847 int nrso = _somap.Extent(), nrsh = _shmap.Extent(), nrf = _fmap.Extent();
4848 int nrw = _wmap.Extent(), nre = _emap.Extent(), nrv = _vmap.Extent();
4849 for(exp0.Init(myshape, TopAbs_COMPOUND); exp0.More(); exp0.Next()) nrc++;
4850 for(exp0.Init(myshape, TopAbs_COMPSOLID); exp0.More(); exp0.Next()) nrcs++;
4851
4852 double surfacecont = 0;
4853 for(exp0.Init(myshape, TopAbs_FACE); exp0.More(); exp0.Next()) {
4854 TopoDS_Face face = TopoDS::Face(exp0.Current());
4855 GProp_GProps system;
4856 BRepGProp::SurfaceProperties(face, system);
4857 surfacecont += system.Mass();
4858 }
4859
4860 if(fixDegenerated) {
4861 Msg::Info(" - Fixing degenerated edges and faces");
4862
4863 {
4864 ShapeBuild_ReShape rebuild;
4865 for(exp1.Init(myshape, TopAbs_EDGE); exp1.More(); exp1.Next()) {
4866 TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
4867 if(BRep_Tool::Degenerated(edge)) rebuild.Remove(edge);
4868 }
4869 myshape = rebuild.Apply(myshape);
4870 }
4871
4872 {
4873 ShapeBuild_ReShape rebuild;
4874 for(exp0.Init(myshape, TopAbs_FACE); exp0.More(); exp0.Next()) {
4875 TopoDS_Face face = TopoDS::Face(exp0.Current());
4876
4877 ShapeFix_Face sff(face);
4878 sff.FixAddNaturalBoundMode() = Standard_True;
4879 sff.FixSmallAreaWireMode() = Standard_True;
4880 sff.Perform();
4881
4882 if(sff.Status(ShapeExtend_DONE1) || sff.Status(ShapeExtend_DONE2) ||
4883 sff.Status(ShapeExtend_DONE3) || sff.Status(ShapeExtend_DONE4) ||
4884 sff.Status(ShapeExtend_DONE5)) {
4885 Msg::Info(" . Repaired face");
4886 if(sff.Status(ShapeExtend_DONE1))
4887 Msg::Info(" . Some wires are fixed");
4888 else if(sff.Status(ShapeExtend_DONE2))
4889 Msg::Info(" . Orientation of wires fixed");
4890 else if(sff.Status(ShapeExtend_DONE3))
4891 Msg::Info(" . Missing seam added");
4892 else if(sff.Status(ShapeExtend_DONE4))
4893 Msg::Info(" . Small area wire removed");
4894 else if(sff.Status(ShapeExtend_DONE5))
4895 Msg::Info(" . Natural bounds added");
4896
4897 TopoDS_Face newface = sff.Face();
4898 rebuild.Replace(face, newface);
4899 }
4900 }
4901 myshape = rebuild.Apply(myshape);
4902 }
4903
4904 {
4905 ShapeBuild_ReShape rebuild;
4906 for(exp1.Init(myshape, TopAbs_EDGE); exp1.More(); exp1.Next()) {
4907 TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
4908 if(BRep_Tool::Degenerated(edge)) rebuild.Remove(edge);
4909 }
4910 myshape = rebuild.Apply(myshape);
4911 }
4912 }
4913
4914 if(fixSmallEdges) {
4915 Msg::Info(" - Fixing small edges");
4916
4917 {
4918 ShapeBuild_ReShape rebuild;
4919
4920 for(exp0.Init(myshape, TopAbs_FACE); exp0.More(); exp0.Next()) {
4921 TopoDS_Face face = TopoDS::Face(exp0.Current());
4922
4923 for(exp1.Init(face, TopAbs_WIRE); exp1.More(); exp1.Next()) {
4924 TopoDS_Wire oldwire = TopoDS::Wire(exp1.Current());
4925 ShapeFix_Wire sfw(oldwire, face, tolerance);
4926 sfw.ModifyTopologyMode() = Standard_True;
4927 sfw.ClosedWireMode() = Standard_True;
4928 bool replace = false;
4929 replace = sfw.FixReorder() || replace;
4930 replace = sfw.FixConnected() || replace;
4931
4932 if(sfw.FixSmall(Standard_False, tolerance) &&
4933 !(sfw.StatusSmall(ShapeExtend_FAIL1) ||
4934 sfw.StatusSmall(ShapeExtend_FAIL2) ||
4935 sfw.StatusSmall(ShapeExtend_FAIL3))) {
4936 Msg::Info(" . Fixed small edge in wire");
4937 replace = true;
4938 }
4939 else if(sfw.StatusSmall(ShapeExtend_FAIL1))
4940 Msg::Warning("Failed to fix small edge in wire, edge cannot be "
4941 "checked (no 3d curve and no pcurve)");
4942 else if(sfw.StatusSmall(ShapeExtend_FAIL2))
4943 Msg::Warning("Failed to fix small edge in wire, edge is null-"
4944 "length and has different vertives at begin and end, "
4945 "and lockvtx is True or ModifiyTopologyMode is False");
4946 else if(sfw.StatusSmall(ShapeExtend_FAIL3))
4947 Msg::Warning("Failed to fix small edge in wire, CheckConnected has "
4948 "failed");
4949
4950 replace = sfw.FixEdgeCurves() || replace;
4951 replace = sfw.FixDegenerated() || replace;
4952 replace = sfw.FixSelfIntersection() || replace;
4953 replace = sfw.FixLacking(Standard_True) || replace;
4954 if(replace) {
4955 TopoDS_Wire newwire = sfw.Wire();
4956 rebuild.Replace(oldwire, newwire);
4957 }
4958 }
4959 }
4960 myshape = rebuild.Apply(myshape);
4961 }
4962
4963 {
4964 ShapeBuild_ReShape rebuild;
4965
4966 for(exp1.Init(myshape, TopAbs_EDGE); exp1.More(); exp1.Next()) {
4967 TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
4968 GProp_GProps system;
4969 BRepGProp::LinearProperties(edge, system);
4970 if(system.Mass() < tolerance) {
4971 Msg::Info(" - Removing degenerated edge");
4972 rebuild.Remove(edge);
4973 }
4974 }
4975 myshape = rebuild.Apply(myshape);
4976 }
4977
4978 {
4979 ShapeBuild_ReShape rebuild;
4980 for(exp1.Init(myshape, TopAbs_EDGE); exp1.More(); exp1.Next()) {
4981 TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
4982 if(BRep_Tool::Degenerated(edge)) rebuild.Remove(edge);
4983 }
4984 myshape = rebuild.Apply(myshape);
4985 }
4986
4987 ShapeFix_Wireframe sfwf;
4988 sfwf.SetPrecision(tolerance);
4989 sfwf.Load(myshape);
4990 sfwf.ModeDropSmallEdges() = Standard_True;
4991
4992 if(sfwf.FixWireGaps()) {
4993 Msg::Info(" - Fixing wire gaps");
4994 if(sfwf.StatusWireGaps(ShapeExtend_OK)) Msg::Info(" no gaps found");
4995 if(sfwf.StatusWireGaps(ShapeExtend_DONE1))
4996 Msg::Info(" . Some 2D gaps fixed");
4997 if(sfwf.StatusWireGaps(ShapeExtend_DONE2))
4998 Msg::Info(" . Some 3D gaps fixed");
4999 if(sfwf.StatusWireGaps(ShapeExtend_FAIL1))
5000 Msg::Info(" . Failed to fix some 2D gaps");
5001 if(sfwf.StatusWireGaps(ShapeExtend_FAIL2))
5002 Msg::Info(" . Failed to fix some 3D gaps");
5003 }
5004
5005 sfwf.SetPrecision(tolerance);
5006
5007 if(sfwf.FixSmallEdges()) {
5008 Msg::Info(" - Fixing wire frames");
5009 if(sfwf.StatusSmallEdges(ShapeExtend_OK))
5010 Msg::Info(" . No small edges found");
5011 if(sfwf.StatusSmallEdges(ShapeExtend_DONE1))
5012 Msg::Info(" . Some small edges fixed");
5013 if(sfwf.StatusSmallEdges(ShapeExtend_FAIL1))
5014 Msg::Info(" . Failed to fix some small edges");
5015 }
5016
5017 myshape = sfwf.Shape();
5018 }
5019
5020 if(fixSmallFaces) {
5021 Msg::Info(" - Fixing spot and strip faces");
5022 ShapeFix_FixSmallFace sffsm;
5023 sffsm.Init(myshape);
5024 sffsm.SetPrecision(tolerance);
5025 sffsm.Perform();
5026 myshape = sffsm.FixShape();
5027 }
5028
5029 if(sewFaces) {
5030 Msg::Info(" - Sewing faces");
5031
5032 BRepOffsetAPI_Sewing sewedObj(tolerance);
5033
5034 for(exp0.Init(myshape, TopAbs_FACE); exp0.More(); exp0.Next()) {
5035 TopoDS_Face face = TopoDS::Face(exp0.Current());
5036 sewedObj.Add(face);
5037 }
5038
5039 sewedObj.Perform();
5040
5041 if(!sewedObj.SewedShape().IsNull())
5042 myshape = sewedObj.SewedShape();
5043 else
5044 Msg::Info(" . Could not sew");
5045 }
5046
5047 {
5048 ShapeBuild_ReShape rebuild;
5049 for(exp1.Init(myshape, TopAbs_EDGE); exp1.More(); exp1.Next()) {
5050 TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
5051 if(BRep_Tool::Degenerated(edge)) rebuild.Remove(edge);
5052 }
5053 myshape = rebuild.Apply(myshape);
5054 }
5055
5056 if(makeSolids) {
5057 Msg::Info(" - Making solids");
5058
5059 BRepBuilderAPI_MakeSolid ms;
5060 int count = 0;
5061 for(exp0.Init(myshape, TopAbs_SHELL); exp0.More(); exp0.Next()) {
5062 count++;
5063 ms.Add(TopoDS::Shell(exp0.Current()));
5064 }
5065
5066 if(!count) { Msg::Info(" . Could not make solid (no shells)"); }
5067 else {
5068 BRepCheck_Analyzer ba(ms);
5069 if(ba.IsValid()) {
5070 ShapeFix_Shape sfs;
5071 sfs.Init(ms);
5072 sfs.SetPrecision(tolerance);
5073 sfs.SetMaxTolerance(tolerance);
5074 sfs.Perform();
5075 myshape = sfs.Shape();
5076 for(exp0.Init(myshape, TopAbs_SOLID); exp0.More(); exp0.Next()) {
5077 TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
5078 TopoDS_Solid newsolid = solid;
5079 BRepLib::OrientClosedSolid(newsolid);
5080 ShapeBuild_ReShape rebuild;
5081 rebuild.Replace(solid, newsolid);
5082 myshape = rebuild.Apply(myshape, TopAbs_COMPSOLID);
5083 }
5084 }
5085 else
5086 Msg::Info(" . Could not make solid");
5087 }
5088 }
5089
5090 double newsurfacecont = 0;
5091 for(exp0.Init(myshape, TopAbs_FACE); exp0.More(); exp0.Next()) {
5092 TopoDS_Face face = TopoDS::Face(exp0.Current());
5093 GProp_GProps system;
5094 BRepGProp::SurfaceProperties(face, system);
5095 newsurfacecont += system.Mass();
5096 }
5097
5098 _somap.Clear();
5099 _shmap.Clear();
5100 _fmap.Clear();
5101 _wmap.Clear();
5102 _emap.Clear();
5103 _vmap.Clear();
5104 _addShapeToMaps(myshape);
5105 int nnrc = 0, nnrcs = 0;
5106 int nnrso = _somap.Extent(), nnrsh = _shmap.Extent(), nnrf = _fmap.Extent();
5107 int nnrw = _wmap.Extent(), nnre = _emap.Extent(), nnrv = _vmap.Extent();
5108 for(exp0.Init(myshape, TopAbs_COMPOUND); exp0.More(); exp0.Next()) nnrc++;
5109 for(exp0.Init(myshape, TopAbs_COMPSOLID); exp0.More(); exp0.Next()) nnrcs++;
5110
5111 double t2 = Cpu(), w2 = TimeOfDay();
5112 Msg::Info("Done healing shapes (Wall %gs, CPU %gs):", w2 - w1, t2 - t1);
5113 Msg::Info(" - Compounds : %d (%d)", nnrc, nrc);
5114 Msg::Info(" - Composite solids : %d (%d)", nnrcs, nrcs);
5115 Msg::Info(" - Solids : %d (%d)", nnrso, nrso);
5116 Msg::Info(" - Shells : %d (%d)", nnrsh, nrsh);
5117 Msg::Info(" - Wires : %d (%d)", nnrw, nrw);
5118 Msg::Info(" - Faces : %d (%d)", nnrf, nrf);
5119 Msg::Info(" - Edges : %d (%d)", nnre, nre);
5120 Msg::Info(" - Vertices : %d (%d)", nnrv, nrv);
5121 Msg::Info(" - Total surface area : %g (%g)", newsurfacecont, surfacecont);
5122 }
5123
healShapes(const std::vector<std::pair<int,int>> & inDimTags,std::vector<std::pair<int,int>> & outDimTags,double tolerance,bool fixDegenerated,bool fixSmallEdges,bool fixSmallFaces,bool sewFaces,bool makeSolids)5124 bool OCC_Internals::healShapes(
5125 const std::vector<std::pair<int, int> > &inDimTags,
5126 std::vector<std::pair<int, int> > &outDimTags, double tolerance,
5127 bool fixDegenerated, bool fixSmallEdges, bool fixSmallFaces, bool sewFaces,
5128 bool makeSolids)
5129 {
5130 BRep_Builder b;
5131 TopoDS_Compound c;
5132 b.MakeCompound(c);
5133
5134 // construct a compound with all the shapes with tags
5135 _somap.Clear();
5136 _shmap.Clear();
5137 _fmap.Clear();
5138 _wmap.Clear();
5139 _emap.Clear();
5140 _vmap.Clear();
5141 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp0(_tagVertex);
5142 for(; exp0.More(); exp0.Next()) _addShapeToMaps(exp0.Value());
5143 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp1(_tagEdge);
5144 for(; exp1.More(); exp1.Next()) _addShapeToMaps(exp1.Value());
5145 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp2(_tagFace);
5146 for(; exp2.More(); exp2.Next()) _addShapeToMaps(exp2.Value());
5147 TopTools_DataMapIteratorOfDataMapOfIntegerShape exp3(_tagSolid);
5148 for(; exp3.More(); exp3.Next()) _addShapeToMaps(exp3.Value());
5149 for(int i = 1; i <= _vmap.Extent(); i++) b.Add(c, _vmap(i));
5150 for(int i = 1; i <= _emap.Extent(); i++) b.Add(c, _emap(i));
5151 for(int i = 1; i <= _wmap.Extent(); i++) b.Add(c, _wmap(i));
5152 for(int i = 1; i <= _fmap.Extent(); i++) b.Add(c, _fmap(i));
5153 for(int i = 1; i <= _shmap.Extent(); i++) b.Add(c, _shmap(i));
5154 for(int i = 1; i <= _somap.Extent(); i++) b.Add(c, _somap(i));
5155
5156 std::vector<TopoDS_Shape> toHeal;
5157 for(std::size_t i = 0; i < inDimTags.size(); i++) {
5158 int dim = inDimTags[i].first;
5159 int tag = inDimTags[i].second;
5160 if(!_isBound(dim, tag)) {
5161 Msg::Error("Unknown OpenCASCADE entity of dimension %d with tag %d",
5162 dim, tag);
5163 return false;
5164 }
5165 toHeal.push_back(_find(dim, tag));
5166 }
5167
5168 // unbind everything
5169 _unbind();
5170
5171 if(toHeal.empty()) { // heal all the shapes
5172 _healShape(c, tolerance, fixDegenerated, fixSmallEdges, fixSmallFaces,
5173 sewFaces, makeSolids, 1.0);
5174 }
5175 else {
5176 for(std::size_t i = 0; i < toHeal.size(); i++) {
5177 TopoDS_Shape olds = toHeal[i];
5178 TopoDS_Shape news = olds;
5179 _healShape(news, tolerance, fixDegenerated, fixSmallEdges, fixSmallFaces,
5180 sewFaces, makeSolids, 1.0);
5181 ShapeBuild_ReShape rebuild;
5182 rebuild.Replace(olds, news);
5183 c = TopoDS::Compound(rebuild.Apply(c));
5184 }
5185 }
5186
5187 // rebind
5188 _multiBind(c, -1, outDimTags, false, true);
5189
5190 return true;
5191 }
5192
makeSTL(const TopoDS_Face & s,std::vector<SPoint2> * verticesUV,std::vector<SPoint3> * verticesXYZ,std::vector<SVector3> * normals,std::vector<int> & triangles)5193 static bool makeSTL(const TopoDS_Face &s, std::vector<SPoint2> *verticesUV,
5194 std::vector<SPoint3> *verticesXYZ,
5195 std::vector<SVector3> *normals, std::vector<int> &triangles)
5196 {
5197 if(CTX::instance()->geom.occDisableSTL) return false;
5198
5199 double lin = CTX::instance()->mesh.stlLinearDeflection;
5200 double ang = CTX::instance()->mesh.stlAngularDeflection;
5201
5202 #if OCC_VERSION_HEX > 0x070300
5203 BRepMesh_IncrementalMesh aMesher(s, lin, Standard_True, ang, Standard_True);
5204 #elif OCC_VERSION_HEX > 0x070000
5205 Bnd_Box aBox;
5206 BRepBndLib::Add(s, aBox);
5207 BRepMesh_FastDiscret::Parameters parameters;
5208 parameters.Deflection = lin;
5209 parameters.Relative = Standard_True;
5210 parameters.Angle = ang;
5211 BRepMesh_FastDiscret aMesher(aBox, parameters);
5212 aMesher.Perform(s);
5213 #else
5214 Bnd_Box aBox;
5215 BRepBndLib::Add(s, aBox);
5216 BRepMesh_FastDiscret aMesher(lin, ang, aBox, Standard_False, Standard_False,
5217 Standard_True, Standard_True);
5218 aMesher.Perform(s);
5219 #endif
5220
5221 TopLoc_Location loc;
5222 Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation(s, loc);
5223
5224 if(triangulation.IsNull()) return false;
5225
5226 if(verticesUV && !triangulation->HasUVNodes()) return false;
5227
5228 if(normals && !triangulation->HasUVNodes()) return false;
5229
5230 int start = 0;
5231 if(verticesUV) start = verticesUV->size();
5232 if(verticesXYZ) start = verticesXYZ->size();
5233 for(int i = 1; i <= triangulation->NbNodes(); i++) {
5234 if(verticesUV) {
5235 #if OCC_VERSION_HEX >= 0x070600
5236 gp_Pnt2d p = triangulation->UVNode(i);
5237 #else
5238 gp_Pnt2d p = (triangulation->UVNodes())(i);
5239 #endif
5240 verticesUV->push_back(SPoint2(p.X(), p.Y()));
5241 }
5242 if(verticesXYZ) {
5243 #if OCC_VERSION_HEX >= 0x070600
5244 gp_Pnt pp = triangulation->Node(i);
5245 #else
5246 gp_Pnt pp = (triangulation->Nodes())(i);
5247 #endif
5248 double x = pp.X(), y = pp.Y(), z = pp.Z();
5249 loc.Transformation().Transforms(x, y, z);
5250 verticesXYZ->push_back(SPoint3(x, y, z));
5251 }
5252 if(normals) {
5253 #if OCC_VERSION_HEX >= 0x070600
5254 gp_Pnt2d p = triangulation->UVNode(i);
5255 #else
5256 gp_Pnt2d p = (triangulation->UVNodes())(i);
5257 #endif
5258 Handle(Geom_Surface) sur = BRep_Tool::Surface(s);
5259 gp_Pnt pnt;
5260 gp_Vec du, dv;
5261 sur->D1(p.X(), p.Y(), pnt, du, dv);
5262 SVector3 t1(du.X(), du.Y(), du.Z());
5263 SVector3 t2(dv.X(), dv.Y(), dv.Z());
5264 SVector3 n(crossprod(t1, t2));
5265 n.normalize();
5266 if(s.Orientation() == TopAbs_REVERSED) n *= -1.;
5267 normals->push_back(n);
5268 }
5269 }
5270 for(int i = 1; i <= triangulation->NbTriangles(); i++) {
5271 #if OCC_VERSION_HEX >= 0x070600
5272 Poly_Triangle triangle = triangulation->Triangle(i);
5273 #else
5274 Poly_Triangle triangle = (triangulation->Triangles())(i);
5275 #endif
5276 int p1, p2, p3;
5277 triangle.Get(p1, p2, p3);
5278 triangles.push_back(start + p1 - 1);
5279 if(s.Orientation() == TopAbs_REVERSED) {
5280 triangles.push_back(start + p3 - 1);
5281 triangles.push_back(start + p2 - 1);
5282 }
5283 else {
5284 triangles.push_back(start + p2 - 1);
5285 triangles.push_back(start + p3 - 1);
5286 }
5287 }
5288 return true;
5289 }
5290
makeEdgeSTLFromFace(const TopoDS_Edge & c,const TopoDS_Face & s,std::vector<SPoint3> * verticesXYZ)5291 bool OCC_Internals::makeEdgeSTLFromFace(const TopoDS_Edge &c,
5292 const TopoDS_Face &s,
5293 std::vector<SPoint3> *verticesXYZ)
5294 {
5295 TopLoc_Location transf;
5296 Handle(Poly_Triangulation) trian = BRep_Tool::Triangulation(s, transf);
5297 if(trian.IsNull()) return false;
5298
5299 Handle(Poly_PolygonOnTriangulation) edgepoly =
5300 BRep_Tool::PolygonOnTriangulation(c, trian, transf);
5301 if(edgepoly.IsNull()) return false;
5302 if(edgepoly->NbNodes() < 2) return false;
5303
5304 for(int i = 1; i <= edgepoly->NbNodes(); i++) {
5305 #if OCC_VERSION_HEX > 0x070600
5306 int j = edgepoly->Node(i);
5307 #else
5308 int j = (edgepoly->Nodes())(i);
5309 #endif
5310 #if OCC_VERSION_HEX >= 0x070600
5311 gp_Pnt pp = trian->Node(j);
5312 #else
5313 gp_Pnt pp = (trian->Nodes())(j);
5314 #endif
5315 if(!transf.IsIdentity()) { pp.Transform(transf); }
5316 verticesXYZ->push_back(SPoint3(pp.X(), pp.Y(), pp.Z()));
5317 }
5318 return true;
5319 }
5320
makeFaceSTL(const TopoDS_Face & s,std::vector<SPoint2> & vertices_uv,std::vector<int> & triangles)5321 bool OCC_Internals::makeFaceSTL(const TopoDS_Face &s,
5322 std::vector<SPoint2> &vertices_uv,
5323 std::vector<int> &triangles)
5324 {
5325 return makeSTL(s, &vertices_uv, nullptr, nullptr, triangles);
5326 }
5327
makeFaceSTL(const TopoDS_Face & s,std::vector<SPoint2> & vertices_uv,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5328 bool OCC_Internals::makeFaceSTL(const TopoDS_Face &s,
5329 std::vector<SPoint2> &vertices_uv,
5330 std::vector<SPoint3> &vertices,
5331 std::vector<SVector3> &normals,
5332 std::vector<int> &triangles)
5333 {
5334 return makeSTL(s, &vertices_uv, &vertices, &normals, triangles);
5335 }
5336
makeFaceSTL(const TopoDS_Face & s,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5337 bool OCC_Internals::makeFaceSTL(const TopoDS_Face &s,
5338 std::vector<SPoint3> &vertices,
5339 std::vector<SVector3> &normals,
5340 std::vector<int> &triangles)
5341 {
5342 return makeSTL(s, nullptr, &vertices, &normals, triangles);
5343 }
5344
_makeSTL(const TopoDS_Shape & s,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5345 bool OCC_Internals::_makeSTL(const TopoDS_Shape &s,
5346 std::vector<SPoint3> &vertices,
5347 std::vector<SVector3> &normals,
5348 std::vector<int> &triangles)
5349 {
5350 bool ret = true;
5351 TopExp_Explorer exp0;
5352 for(exp0.Init(s, TopAbs_FACE); exp0.More(); exp0.Next()) {
5353 TopoDS_Face face = TopoDS::Face(exp0.Current());
5354 bool tmp = makeSTL(TopoDS::Face(face.Oriented(TopAbs_FORWARD)), nullptr,
5355 &vertices, &normals, triangles);
5356 if(!tmp) ret = false;
5357 }
5358 return ret;
5359 }
5360
makeRectangleSTL(double x,double y,double z,double dx,double dy,double roundedRadius,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5361 bool OCC_Internals::makeRectangleSTL(double x, double y, double z, double dx,
5362 double dy, double roundedRadius,
5363 std::vector<SPoint3> &vertices,
5364 std::vector<SVector3> &normals,
5365 std::vector<int> &triangles)
5366 {
5367 TopoDS_Face result;
5368 if(!makeRectangle(result, x, y, z, dx, dy, roundedRadius)) return false;
5369 if(!makeFaceSTL(result, vertices, normals, triangles)) return false;
5370 return true;
5371 }
5372
makeDiskSTL(double xc,double yc,double zc,double rx,double ry,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5373 bool OCC_Internals::makeDiskSTL(double xc, double yc, double zc, double rx,
5374 double ry, std::vector<SPoint3> &vertices,
5375 std::vector<SVector3> &normals,
5376 std::vector<int> &triangles)
5377 {
5378 TopoDS_Face result;
5379 if(!makeDisk(result, xc, yc, zc, rx, ry)) return false;
5380 if(!makeFaceSTL(result, vertices, normals, triangles)) return false;
5381 return true;
5382 }
5383
makeSphereSTL(double xc,double yc,double zc,double radius,double angle1,double angle2,double angle3,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5384 bool OCC_Internals::makeSphereSTL(double xc, double yc, double zc,
5385 double radius, double angle1, double angle2,
5386 double angle3, std::vector<SPoint3> &vertices,
5387 std::vector<SVector3> &normals,
5388 std::vector<int> &triangles)
5389 {
5390 TopoDS_Solid result;
5391 if(!makeSphere(result, xc, yc, zc, radius, angle1, angle2, angle3))
5392 return false;
5393 if(!_makeSTL(result, vertices, normals, triangles)) return false;
5394 return true;
5395 }
5396
makeBoxSTL(double x,double y,double z,double dx,double dy,double dz,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5397 bool OCC_Internals::makeBoxSTL(double x, double y, double z, double dx,
5398 double dy, double dz,
5399 std::vector<SPoint3> &vertices,
5400 std::vector<SVector3> &normals,
5401 std::vector<int> &triangles)
5402 {
5403 TopoDS_Solid result;
5404 if(!makeBox(result, x, y, z, dx, dy, dz)) return false;
5405 if(!_makeSTL(result, vertices, normals, triangles)) return false;
5406 return true;
5407 }
5408
makeCylinderSTL(double x,double y,double z,double dx,double dy,double dz,double r,double angle,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5409 bool OCC_Internals::makeCylinderSTL(double x, double y, double z, double dx,
5410 double dy, double dz, double r,
5411 double angle,
5412 std::vector<SPoint3> &vertices,
5413 std::vector<SVector3> &normals,
5414 std::vector<int> &triangles)
5415 {
5416 TopoDS_Solid result;
5417 if(!makeCylinder(result, x, y, z, dx, dy, dz, r, angle)) return false;
5418 if(!_makeSTL(result, vertices, normals, triangles)) return false;
5419 return true;
5420 }
5421
makeConeSTL(double x,double y,double z,double dx,double dy,double dz,double r1,double r2,double angle,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5422 bool OCC_Internals::makeConeSTL(double x, double y, double z, double dx,
5423 double dy, double dz, double r1, double r2,
5424 double angle, std::vector<SPoint3> &vertices,
5425 std::vector<SVector3> &normals,
5426 std::vector<int> &triangles)
5427 {
5428 TopoDS_Solid result;
5429 if(!makeCone(result, x, y, z, dx, dy, dz, r1, r2, angle)) return false;
5430 if(!_makeSTL(result, vertices, normals, triangles)) return false;
5431 return true;
5432 }
5433
makeWedgeSTL(double x,double y,double z,double dx,double dy,double dz,double ltx,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5434 bool OCC_Internals::makeWedgeSTL(double x, double y, double z, double dx,
5435 double dy, double dz, double ltx,
5436 std::vector<SPoint3> &vertices,
5437 std::vector<SVector3> &normals,
5438 std::vector<int> &triangles)
5439 {
5440 TopoDS_Solid result;
5441 if(!makeWedge(result, x, y, z, dx, dy, dz, ltx)) return false;
5442 if(!_makeSTL(result, vertices, normals, triangles)) return false;
5443 return true;
5444 }
5445
makeTorusSTL(double x,double y,double z,double r1,double r2,double angle,std::vector<SPoint3> & vertices,std::vector<SVector3> & normals,std::vector<int> & triangles)5446 bool OCC_Internals::makeTorusSTL(double x, double y, double z, double r1,
5447 double r2, double angle,
5448 std::vector<SPoint3> &vertices,
5449 std::vector<SVector3> &normals,
5450 std::vector<int> &triangles)
5451 {
5452 TopoDS_Solid result;
5453 if(!makeTorus(result, x, y, z, r1, r2, angle)) return false;
5454 if(!_makeSTL(result, vertices, normals, triangles)) return false;
5455 return true;
5456 }
5457
fixSTLBounds(double & xmin,double & ymin,double & zmin,double & xmax,double & ymax,double & zmax)5458 void OCC_Internals::fixSTLBounds(double &xmin, double &ymin, double &zmin,
5459 double &xmax, double &ymax, double &zmax)
5460 {
5461 // When an STL exists, OCC enlarges the bounding box by the allowed linear
5462 // deflection given to BRepMesh_IncrementalMesh. This is "safe", but on simple
5463 // polyhedral geometries (a cube!) it will consistently lead to enlarging the
5464 // bounding box by twice this value in all directions. Since we use bounds()
5465 // mostly for locating entities, it's better to remove the tolerance (with the
5466 // risk that the bbox is a bit too small for curved boundaries - but that's
5467 // fine)
5468 double eps = CTX::instance()->mesh.stlLinearDeflection;
5469 // OCC also enlarges the bounding box by Precision::Confusion(): remove it as
5470 // well
5471 eps += Precision::Confusion();
5472 xmin += eps;
5473 xmax -= eps;
5474 ymin += eps;
5475 ymax -= eps;
5476 zmin += eps;
5477 zmax -= eps;
5478 }
5479
5480 #endif
5481
createOCCInternals()5482 void GModel::createOCCInternals()
5483 {
5484 if(_occ_internals) delete _occ_internals;
5485 _occ_internals = new OCC_Internals;
5486 }
5487
deleteOCCInternals()5488 void GModel::deleteOCCInternals()
5489 {
5490 if(_occ_internals) delete _occ_internals;
5491 _occ_internals = nullptr;
5492 }
5493
resetOCCInternals()5494 void GModel::resetOCCInternals()
5495 {
5496 if(!_occ_internals) return;
5497 _occ_internals->reset();
5498 }
5499
readOCCBREP(const std::string & fn)5500 int GModel::readOCCBREP(const std::string &fn)
5501 {
5502 if(!_occ_internals) _occ_internals = new OCC_Internals;
5503 std::vector<std::pair<int, int> > outDimTags;
5504 _occ_internals->importShapes(fn, false, outDimTags, "brep");
5505 _occ_internals->synchronize(this);
5506 snapVertices();
5507 return 1;
5508 }
5509
readOCCSTEP(const std::string & fn)5510 int GModel::readOCCSTEP(const std::string &fn)
5511 {
5512 if(!_occ_internals) _occ_internals = new OCC_Internals;
5513 std::vector<std::pair<int, int> > outDimTags;
5514 _occ_internals->importShapes(fn, false, outDimTags, "step");
5515 _occ_internals->synchronize(this);
5516 return 1;
5517 }
5518
readOCCIGES(const std::string & fn)5519 int GModel::readOCCIGES(const std::string &fn)
5520 {
5521 if(!_occ_internals) _occ_internals = new OCC_Internals;
5522 std::vector<std::pair<int, int> > outDimTags;
5523 _occ_internals->importShapes(fn, false, outDimTags, "iges");
5524 _occ_internals->synchronize(this);
5525 return 1;
5526 }
5527
writeOCCBREP(const std::string & fn)5528 int GModel::writeOCCBREP(const std::string &fn)
5529 {
5530 if(!_occ_internals) {
5531 Msg::Error("No OpenCASCADE model found");
5532 return 0;
5533 }
5534 _occ_internals->exportShapes(this, fn, "brep");
5535 return 1;
5536 }
5537
writeOCCSTEP(const std::string & fn)5538 int GModel::writeOCCSTEP(const std::string &fn)
5539 {
5540 if(!_occ_internals) {
5541 Msg::Error("No OpenCASCADE model found");
5542 return 0;
5543 }
5544 _occ_internals->exportShapes(this, fn, "step");
5545 return 1;
5546 }
5547
importOCCShape(const void * shape)5548 int GModel::importOCCShape(const void *shape)
5549 {
5550 if(!_occ_internals) _occ_internals = new OCC_Internals;
5551 #if defined(HAVE_OCC)
5552 std::vector<std::pair<int, int> > outDimTags;
5553 _occ_internals->importShapes((TopoDS_Shape *)shape, false, outDimTags);
5554 #else
5555 Msg::Error("Gmsh requires OpenCASCADE to import TopoDS_Shape");
5556 #endif
5557 _occ_internals->synchronize(this);
5558 snapVertices();
5559 return 1;
5560 }
5561
getVertexForOCCShape(const void * shape)5562 GVertex *GModel::getVertexForOCCShape(const void *shape)
5563 {
5564 if(!_occ_internals) return nullptr;
5565 #if defined(HAVE_OCC)
5566 return _occ_internals->getVertexForOCCShape(this, *(TopoDS_Vertex *)shape);
5567 #else
5568 return 0;
5569 #endif
5570 }
5571
getEdgeForOCCShape(const void * shape)5572 GEdge *GModel::getEdgeForOCCShape(const void *shape)
5573 {
5574 if(!_occ_internals) return nullptr;
5575 #if defined(HAVE_OCC)
5576 return _occ_internals->getEdgeForOCCShape(this, *(TopoDS_Edge *)shape);
5577 #else
5578 return 0;
5579 #endif
5580 }
5581
getFaceForOCCShape(const void * shape)5582 GFace *GModel::getFaceForOCCShape(const void *shape)
5583 {
5584 if(!_occ_internals) return nullptr;
5585 #if defined(HAVE_OCC)
5586 return _occ_internals->getFaceForOCCShape(this, *(TopoDS_Face *)shape);
5587 #else
5588 return 0;
5589 #endif
5590 }
5591
getRegionForOCCShape(const void * shape)5592 GRegion *GModel::getRegionForOCCShape(const void *shape)
5593 {
5594 if(!_occ_internals) return nullptr;
5595 #if defined(HAVE_OCC)
5596 return _occ_internals->getRegionForOCCShape(this, *(TopoDS_Solid *)shape);
5597 #else
5598 return 0;
5599 #endif
5600 }
5601