G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsWriteParasolid.hpp
Go to the documentation of this file.
1
14#include <gsParasolid/gsFrustrum.h>
15
16#include <gsParasolid/gsReadParasolid.h>
17#include <gsParasolid/gsPKSession.h>
18
19#include <gsCore/gsMultiPatch.h>
20#include <gsCore/gsSurface.h>
21
22#include <gsNurbs/gsKnotVector.h>
23#include <gsNurbs/gsBSplineBasis.h>
24#include <gsNurbs/gsTensorBSpline.h>
25#include <gsNurbs/gsBSpline.h>
26
27#include <gsUtils/gsMesh/gsMesh.h>
28
29#include <gsHSplines/gsTHBSpline.h>
30#include <gsHSplines/gsTHBSplineBasis.h>
31
32#include <gsIO/gsWriteParaview.h>
33
34
35namespace gismo {
36
37namespace extensions {
38
39/*
40 Notes:
41
42// The sheet is a kind of part
43// PK_PART_t part = sheet;
44
45This function attaches surfaces to faces.
46PK_FACE_attach_surfs
47
48//This function creates a sheet body from a collection of faces.
49PK_FACE_make_sheet_body
50PK_FACE_make_sheet_bodies
51
52//This function creates solid bodies from a collection of faces
53PK_FACE_make_solid_bodies
54
55//PK_FACE_make_sheet_body
56//PK_FACE_make_solid_bodies
57
58//PK_PART_find_entity_by_ident
59
60
61// Create parasolid sheet-body out of the surface
62//PK_UVBOX_t uvbox; // Parameter domain intevals in u and v
63//PK_SURF_ask_uvbox(bsurf, &uvbox); // Get the parameter box
64//PK_BODY_t sheet; // Sheet is a kind of body
65//err = PK_SURF_make_sheet_body(bsurf, uvbox, &sheet);
66//PARASOLID_ERROR(PK_SURF_make_sheet_body, err);
67
68*/
69
70template<class T>
71class gsTrimData
72{
73public:
74 gsTrimData(unsigned level,std::vector<index_t>& AABBBox,std::vector<std::vector<std::vector<T> > >polylines)
75 : m_level(level),m_AABBBox(AABBBox),m_Polylines(polylines)
76 { }
77
78 gsTrimData() : m_level(0),m_AABBBox(),m_Polylines()
79 { }
80
81public:
82 unsigned m_level;
83 std::vector<index_t> m_AABBBox;
84 std::vector<std::vector<std::vector<T> > > m_Polylines;
85};
86
87// forward declarations
88void makeValidGeometry(const gsTHBSpline<2>& surface,
89 gsTensorBSpline<2, real_t>& bspline);
90
91void getInterval(const bool directionU,
92 const real_t param1,
93 const real_t param2,
94 const real_t paramConst,
95 const gsTensorBSpline<2, real_t>& bspline,
96 const PK_CURVE_t line,
97 PK_INTERVAL_t& result);
98
99bool validMultiplicities(const std::vector<index_t>& mult,
100 const int deg);
101
102template <class T>
103bool exportCheck(const gsTHBSpline<2, T>& surface);
104
105template <class T>
106bool exportTHBsurface( const gsTHBSpline<2, T>& surface,
107 gsTHBSplineBasis<2>::TrimmingCurves trimCurves,
108 gsTHBSplineBasis<2>::AxisAlignedBoundingBox boundaryAABB,
109 PK_ASSEMBLY_t& body );
110
111template <class T>
112void getTrimCurvesAndBoundingBoxes(const gsTHBSpline<2, T>& surface,
113 std::vector<index_t>& boxes,
114 gsTHBSplineBasis<2>::TrimmingCurves& trimCurves,
115 gsTHBSplineBasis<2>::AxisAlignedBoundingBox& boundaryAABB);
116
117template<class T>
118bool gsWriteParasolid( const gsMultiPatch<T> & gssurfs, std::string const & filename )
119{
120 std::cout << "write parasolid mulitpatch" << std::endl;
121
122 PK_ERROR_code_t err;
123 PK_BODY_t part; // Empty part
124 PK_GEOM_t geo[ gssurfs.nPatches() ]; // Geometries
125
126 // Start Parasolid session
127 gsPKSession::start();
128
129 // Disable continuity and self-intersection checks on geometrical
130 // data of bodies
131 PK_LOGICAL_t checks(0);
132 PK_SESSION_set_check_continuity(checks);
133 PK_SESSION_set_check_self_int(checks);
134
135 // Create Parasolid geometries
136 int count = 0;
137 for (typename gsMultiPatch<T>::const_iterator
138 it = gssurfs.begin(); it != gssurfs.end(); ++it)
139 {
140 createPK_GEOM(**it, geo[count++] );
141 }
142
143/*
144// Push all geometries as orphans in an assembly part
145PK_ASSEMBLY_create_empty(&part);
146err = PK_PART_add_geoms(part, count, geo);
147PARASOLID_ERROR(PK_PART_add_geoms, err);
148*/
149
150 // Make a sheet body out of each geometry
151 PK_BODY_t parts[count];
152 for (int i = 0; i!= count; ++i )
153 {
154 PK_UVBOX_t uvbox; // Parameter interval
155 PK_ERROR_code_t err = PK_SURF_ask_uvbox(geo[i], &uvbox);
156 PARASOLID_ERROR(PK_SURF_ask_uvbox, err);
157
158 PK_SURF_make_sheet_body(geo[i], uvbox, &part);
159 PARASOLID_ERROR(PK_SURF_make_sheet_body, err);
160
161 parts[i] = part;
162 }
163
164 // Write it out to the file
165 PK_PART_transmit_o_t transmit_options;
166 PK_PART_transmit_o_m(transmit_options);
167 transmit_options.transmit_format = PK_transmit_format_text_c;
168 err = PK_PART_transmit(count, parts, filename.c_str(), &transmit_options);
169 PARASOLID_ERROR(PK_PART_transmit, err);
170
171 // Stop Parasolid session
172 gsPKSession::stop();
173
174 return err;
175}
176
177template<class T>
178bool gsWriteParasolid( const gsGeometry<T> & ggeo, std::string const & filename )
179{
180 PK_ERROR_code_t err;
181 PK_ASSEMBLY_t part;
182 PK_GEOM_t pgeo; // Parasolid geometric entity
183
184 // Start Parasolid session
185 gsPKSession::start();
186
187 // Disable continuity and self-intersectiion checks
188 PK_LOGICAL_t checks(0);
189 PK_SESSION_set_check_continuity(checks);
190 PK_SESSION_set_check_self_int(checks);
191
192 // Create Parasolid geometry
193 createPK_GEOM(ggeo, pgeo);
194
195 // Create parasolid part out of the geometry
196 PK_ASSEMBLY_create_empty(&part);
197 err = PK_PART_add_geoms(part, 1, &pgeo);
198 PARASOLID_ERROR(PK_PART_add_geoms, err);
199
200 // Write it out to the file
201 PK_PART_transmit_o_t transmit_options;
202 PK_PART_transmit_o_m(transmit_options);
203 transmit_options.transmit_format = PK_transmit_format_text_c;
204 err = PK_PART_transmit(1, &part, filename.c_str(), &transmit_options);
205 PARASOLID_ERROR(PK_PART_transmit, err);
206
207 // Stop Parasolid session
208 gsPKSession::stop();
209
210 return err;
211}
212
213
214template<class T>
215bool gsWriteParasolid( const gsMesh<T>& mesh, const std::string & filename)
216{
217 gsPKSession::start();
218
219 PK_LOGICAL_t checks(0);
220 PK_SESSION_set_check_continuity(checks);
221 PK_SESSION_set_check_self_int(checks);
222
223 PK_ASSEMBLY_t assembly;
224 exportMesh(mesh, assembly);
225
226 PK_PART_transmit_o_t transmit_options;
227 PK_PART_transmit_o_m(transmit_options);
228 transmit_options.transmit_format = PK_transmit_format_text_c;
229
230 PK_ERROR_code_t err = PK_PART_transmit(1, &assembly, filename.c_str(), &transmit_options);
231 PARASOLID_ERROR(PK_PART_transmit, err);
232
233 gsPKSession::stop();
234
235 return err;
236}
237
238
239template<class T>
240bool gsWriteParasolid(const gsTHBSpline<2, T>& thb, const std::string& filename)
241{
242 gsPKSession::start();
243 PK_LOGICAL_t checks(0);
244 PK_SESSION_set_check_continuity(checks);
245 PK_SESSION_set_check_self_int(checks);
246
247 PK_ASSEMBLY_t assembly;
248 exportTHBsurface<T>(thb, assembly);
249
250 PK_PART_transmit_o_t transmit_options;
251 PK_PART_transmit_o_m(transmit_options);
252 transmit_options.transmit_format = PK_transmit_format_text_c;
253
254 PK_ERROR_code_t err = PK_PART_transmit(1, &assembly, filename.c_str(), &transmit_options);
255 PARASOLID_ERROR(PK_PART_transmit, err);
256
257 gsPKSession::stop();
258
259 return err;
260}
261
262template<class T>
263bool gsWriteParasolid( const gsTHBSpline<2, T>& thb,const std::vector<T>& par_boxes, const std::string& filename )
264{
265 gsPKSession::start();
266 PK_LOGICAL_t checks(0);
267 PK_SESSION_set_check_continuity(checks);
268 PK_SESSION_set_check_self_int(checks);
269
270 PK_ASSEMBLY_t assembly;
271 exportTHBsurface<T>(thb,par_boxes,assembly);
272
273 PK_PART_transmit_o_t transmit_options;
274 PK_PART_transmit_o_m(transmit_options);
275 transmit_options.transmit_format = PK_transmit_format_text_c;
276
277 PK_ERROR_code_t err = PK_PART_transmit(1, &assembly, filename.c_str(), &transmit_options);
278 PARASOLID_ERROR(PK_PART_transmit, err);
279 gsPKSession::stop();
280 return err;
281}
282
283template<class T>
284bool gsWritePK_SHEET(const gsTensorBSpline<2, T>& tp, const std::string& filename)
285{
286 gsPKSession::start();
287 PK_LOGICAL_t checks(0);
288 PK_SESSION_set_check_continuity(checks);
289 PK_SESSION_set_check_self_int(checks);
290
291 PK_ERROR_code_t err;
292
293 PK_BSURF_t bsurf;
294 createPK_BSURF<T>(tp, bsurf, false, false);
295
296 gsMatrix<T> domain = tp.parameterRange();
297 PK_UVBOX_t uv_box;
298 // format: umin, vmin, umax, vmax
299 uv_box.param[0] = domain(0, 0);
300 uv_box.param[1] = domain(1, 0);
301 uv_box.param[2] = domain(0, 1);
302 uv_box.param[3] = domain(1, 1);
303 PK_BODY_t body;
304 err = PK_SURF_make_sheet_body(bsurf, uv_box, &body);
305 PARASOLID_ERROR(PK_SURF_make_sheet_body, err);
306
307 // PK_ASSEMBLY_t assembly;
308 // PK_ASSEMBLY_create_empty(&assembly);
309 // err = PK_PART_add_geoms(assembly, 1, &bsurf);
310 // PARASOLID_ERROR(PK_PART_add_geoms, err);
311
312 PK_PART_transmit_o_t transmit_options;
313 PK_PART_transmit_o_m(transmit_options);
314 transmit_options.transmit_format = PK_transmit_format_text_c;
315
316 err = PK_PART_transmit(1, &body, filename.c_str(), &transmit_options);
317 PARASOLID_ERROR(PK_PART_transmit, err);
318
319 gsPKSession::stop();
320
321 return err;
322}
323
324
325template<class T>
326bool createPK_GEOM( const gsGeometry<T> & ggeo,
327 PK_GEOM_t & pgeo)
328{
329 // Identify input gismo geometry
330 if ( const gsTensorBSpline<2,T> * tbsp =
331 dynamic_cast<const gsTensorBSpline<2,T> *>(&ggeo) )
332 {
333 return createPK_BSURF(*tbsp, pgeo);
334 }
335// the following lines produce warnings if called from multipatch version of gsWriteParasolid,
336// because it already assumes that the geometries are surfaces
337 else if ( const gsBSpline<>* bspl =
338 dynamic_cast< const gsBSpline<>* >(&ggeo) )
339 {
340 return createPK_BCURVE(*bspl, pgeo);
341 }
342 else
343 {
344 gsInfo << "Cannot write "<<ggeo<<" to parasolid file.\n";
345 return false;
346 }
347}
348
349
350template<class T>
351bool createPK_BSURF(const gsTensorBSpline< 2, T> & bsp,
352 PK_BSURF_t & bsurf,
353 bool closed_u,
354 bool closed_v)
355{
356 typedef typename gsKnotVector<T>::mult_t mult_t;
357 std::vector<mult_t> mult;
358
359 // Gismo and Parasolid store the coefficients in different order.
360 // Therefore we create the BSURF with u and v swapped and then transpose it.
361 for (index_t dim = 0; dim != 2; dim++)
362 {
363 const int deg = bsp.basis().degree(dim);
364 mult = bsp.basis().knots(dim).multiplicities();
365
366 if (!validMultiplicities(mult, deg))
367 {
368 return false;
369 }
370 }
371
372 // Translate to parasolid standard form, ie fill up parasolid
373 // spline data record
374 PK_BSURF_sf_t sform; // B-spline data holder (standard form)
375
376 // Degrees
377 sform.u_degree = bsp.basis().degree(1);
378 sform.v_degree = bsp.basis().degree(0);
379
380 // Knots in u-direction
381 std::vector<T> gknot0 = bsp.basis().knots(1).unique();
382 std::vector<mult_t> gmult0 = bsp.basis().knots(1).multiplicities();
383 sform.n_u_knots = gknot0.size();
384 sform.u_knot = gknot0.data();
385 sform.u_knot_mult = gmult0.data();
386
387 // Knots in v-direction
388 std::vector<T> gknot1 = bsp.basis().knots(0).unique();
389 std::vector<mult_t> gmult1 = bsp.basis().knots(0).multiplicities();
390 sform.n_v_knots = gknot1.size();
391 sform.v_knot = gknot1.data();
392 sform.v_knot_mult = gmult1.data();
393
394 // Control points
395 sform.n_u_vertices = bsp.basis().size(1);
396 sform.n_v_vertices = bsp.basis().size(0);
397 gsMatrix<T> coefs = bsp.coefs();
398 const int n = bsp.geoDim();
399 if ( n < 3 )
400 {
401 coefs.conservativeResize(gsEigen::NoChange, 3);
402 coefs.rightCols(3-n).setZero();
403 }
404 coefs.transposeInPlace();
405 coefs.resize(3*bsp.basis().size(), 1);
406 sform.vertex_dim = 3; // always 3 for surfaces
407 sform.vertex = coefs.data();
408
409 // Attributes
410 sform.is_rational = PK_LOGICAL_false;
411 sform.form = PK_BSURF_form_unset_c;
412 sform.u_knot_type = PK_knot_unset_c;
413 sform.v_knot_type = PK_knot_unset_c;
414 sform.is_u_periodic = PK_LOGICAL_false;
415 sform.is_v_periodic = PK_LOGICAL_false;
416 sform.is_u_closed = PK_LOGICAL_false;
417 sform.is_v_closed = PK_LOGICAL_false;
418
419 if (closed_u)
420 {
421 sform.is_v_closed = PK_LOGICAL_true;
422 }
423
424 if (closed_v)
425 {
426 sform.is_u_closed = PK_LOGICAL_true;
427 }
428
429 sform.self_intersecting = PK_self_intersect_unset_c;
430 sform.convexity = PK_convexity_unset_c;
431
432 // Create parasolid surface with the previous spline data
433 PK_ERROR_code_t err = PK_BSURF_create(&sform, &bsurf);
434 PARASOLID_ERROR(PK_BSURF_create, err);
435
436 // Transposition (new on 2019-02-26).
437 PK_BSURF_reparameterise_o_t options;
438 PK_BSURF_reparameterise_o_m(options);
439 options.transpose = PK_LOGICAL_true;
440
441 err = PK_BSURF_reparameterise(bsurf,&options);
442 PARASOLID_ERROR(PK_BSURF_reparameterise, err);
443
444 return true;
445}
446
447template<class T>
448bool createPK_BCURVE( const gsBSpline<T>& curve,
449 PK_BCURVE_t& bcurve)
450{
451 typedef typename gsKnotVector<T>::mult_t mult_t;
452
453 PK_BCURVE_sf_t sform; // B-curve data holder (standard form)
454
455 // Degree
456 sform.degree = curve.degree();
457
458 // Knots
459 std::vector<T> knots = curve.basis().knots().unique();
460 std::vector<mult_t> mult = curve.basis().knots().multiplicities();
461 sform.n_knots = knots.size();
462 sform.knot = knots.data();
463 sform.knot_mult = mult.data();
464
465
466 // Control points
467 sform.n_vertices = curve.basis().size();
468 gsMatrix<T> coefs = curve.coefs();
469 const int n = curve.geoDim();
470 if (n < 3)
471 {
472 coefs.conservativeResize(gsEigen::NoChange, 3);
473 coefs.rightCols(3 - n).setZero();
474 }
475 coefs.transposeInPlace();
476 coefs.resize(3 * curve.basis().size(), 1);
477 sform.vertex_dim = 3;
478 sform.vertex = coefs.data();
479
480 // Attributes
481 sform.is_rational = PK_LOGICAL_false;
482 sform.form = PK_BCURVE_form_unset_c;
483 sform.knot_type = PK_knot_unset_c;
484 sform.is_periodic = PK_LOGICAL_false;
485 sform.is_closed = PK_LOGICAL_false;
486 sform.self_intersecting = PK_self_intersect_unset_c;
487
488 PK_ERROR_code_t err = PK_BCURVE_create(&sform, &bcurve);
489 PARASOLID_ERROR(PK_BCURVE_create, err);
490
491 return true;
492}
493
494template<class T>
495bool exportMesh(const gsMesh<T>& mesh,
496 PK_ASSEMBLY_t& assembly)
497{
498 // tried to:
499 // - make one wire body out of all mesh edges with PK_CURVE_make_wire_body_2
500 // and it doesn't work, because edges cross each other
501 // - make a boolean union of all mesh edges with PK_BODY_boolean_2
502 // and it doesn't work function fails to make a union body (I don't know
503 // the reason)
504
505 PK_ERROR_t err = PK_ASSEMBLY_create_empty(&assembly);
506 PARASOLID_ERROR(PK_ASSEMBLY_create_empty, err);
507
508 gsKnotVector<T> kv(0, 1, 0, 2);
509 gsMatrix<T> coefs(2, 3);
510 gsBSpline<T> bspl(kv, coefs);
511
512 gsMatrix<T> newCoefs(2, 3);
513 for (size_t i = 0; i != mesh.numEdges(); ++i)
514 {
515 newCoefs.row(0) = mesh.edges()[i].source->transpose();
516 newCoefs.row(1) = mesh.edges()[i].target->transpose();
517
518 if ((newCoefs.row(0) - newCoefs.row(1)).norm() < 1e-6)
519 {
520 continue;
521 }
522
523 bspl.setCoefs(newCoefs);
524
525 PK_BCURVE_t bcurve;
526 createPK_BCURVE(bspl, bcurve);
527
528 err = PK_PART_add_geoms(assembly, 1, &bcurve);
529 PARASOLID_ERROR(PK_PART_add_geoms, err);
530 }
531
532 return true;
533}
534
535
536template <class T>
537bool exportTHBsurface(const gsTHBSpline<2, T>& surface,
538 PK_ASSEMBLY_t& assembly)
539{
540 if(!exportCheck(surface))
541 return false;
542
543 gsTHBSplineBasis<2>::AxisAlignedBoundingBox boundaryAABB;
544 gsTHBSplineBasis<2>::TrimmingCurves trimCurves;
545
546 const gsTHBSplineBasis<2>& basis= surface.basis();
547
548 basis.decomposeDomain(boundaryAABB, trimCurves);
549
550 std::vector<gsTrimData<T> > trimData;
551
552 bool success = getTrimDataFromBoundaryTrimCurves(boundaryAABB,trimCurves,trimData);
553
554 if(!success)
555 return false;
556
557 return exportTHBsurface(surface,trimData,assembly);
558}
559
560template <class T>
561bool exportTHBsurface(const gsTHBSpline<2, T>& surface,
562 const std::vector<T>& par_boxes,
563 PK_ASSEMBLY_t& assembly)
564{
565 if(par_boxes.size()==0)
566 return exportTHBsurface(surface,assembly);
567
568 if(!exportCheck(surface))
569 return false;
570
571 std::vector<gsTrimData<T> > trimData;
572
573 bool success = getTrimCurvesAndBoundingBoxes<T>(surface,par_boxes,trimData);
574
575 if(!success)
576 return false;
577
578 return exportTHBsurface(surface,trimData,assembly);
579}
580
581
582// ================================================================================
583// utilities
584
585// returns true la all multiplicities in muls are less (<) than deg + 1
586// parasolid restriction
587bool validMultiplicities(const std::vector<index_t>& mult,
588 const int deg)
589{
590 for (size_t i = 1; i != mult.size() - 1; i++)
591 {
592 if (mult[i] == deg + 1)
593 {
594 std::cout <<
595 "Only B-Splines with (inner) multiplicity less (<) "
596 "than degree + 1 are supported. \n"
597 "Parasolid restriction."
598 << std::endl;
599 return false;
600 }
601 }
602 return true;
603}
604
605
606// computes the parameter interval of the iso-curve (line) on the surface (bspline)
607// curve is between two points on the surface defined by input parameters
608// - directionU tells if the curve is is-curve in parametric directionU
609// - paramConst if the fixed parameter in parametric direction defined by directionU
610// - param1 & param2 are non fixed parameters defined in the other direction than before
611// - param1, param2, paramConst defines two points -- between these two points we would like
612// to define parameter range for our curve
613void getInterval(const bool directionU,
614 const real_t param1,
615 const real_t param2,
616 const real_t paramConst,
617 const gsTensorBSpline<2, real_t>& bspline,
618 const PK_CURVE_t line,
619 PK_INTERVAL_t& result)
620{
621 gsMatrix<> param(2, 2);
622 if (directionU)
623 {
624 param(0, 0) = param1;
625 param(1, 0) = paramConst;
626
627 param(0, 1) = param2;
628 param(1, 1) = paramConst;
629 }
630 else
631 {
632 param(0, 0) = paramConst;
633 param(1, 0) = param1;
634
635 param(0, 1) = paramConst;
636 param(1, 1) = param2;
637 }
638
639 gsMatrix<> points;
640 bspline.eval_into(param, points);
641
642
643 for (index_t i = 0; i != 2; i++)
644 {
645 PK_VECTOR_t position;
646 position.coord[0] = 0;
647 position.coord[1] = 0;
648 position.coord[2] = 0;
649
650 for (index_t row = 0; row != points.rows(); row++)
651 {
652 position.coord[row] = points(row, i);
653 }
654
655 real_t p;
656
657 PK_ERROR_t err = PK_CURVE_parameterise_vector(line, position, &p);
658 PARASOLID_ERROR(PK_CURVE_parameterise_vector, err);
659
660 result.value[i] = p;
661 }
662
663 // try to correct the problem of periodic surfaces
664 if( param1!=param2 && result.value[0]==result.value[1] )
665 {
666 result.value[0]=param1;
667 result.value[1]=param2;
668 }
669}
670
671
672// This function is here as workaround around Parasolid / MTU visualization bug.
673// The bug is:
674// Let say that you make a trimmed sheet in a such way that you make a
675// hole into geometry. So there is a region in geometry which will be trimmed
676// away. If multiple coefficients, which define this region, are equal to
677// (0, 0, 0) then Parasolid / MTU visualization complains because the surface
678// self intersects. Parasolid / MTU visulation doesn't realize that this reigon
679// will be trimmed away and it is not important.
680//
681// The workaround:
682// We change the problematic coefficients, from the trimmed-away area, such that
683// the bspline approximates the THB surface.
684//
685// TODO: refactor this function: split it to multiple parts
686void
687makeValidGeometry(const gsTHBSpline<2>& surface,
688 gsTensorBSpline<2, real_t>& bspline)
689{
690 // check how many coefs are zero
691
692 gsMatrix<>& coefs = bspline.coefs();
693 gsVector<int> globalToLocal(coefs.rows());
694 globalToLocal.setConstant(-1);
695
696 int numZeroRows = 0;
697 for (index_t row = 0; row != coefs.rows(); row++)
698 {
699 bool zeroRow = true;
700 for (index_t col = 0; col != coefs.cols(); col++)
701 {
702 if (coefs(row, col) != 0.0)
703 {
704 zeroRow = false;
705 break;
706 }
707 }
708
709 if (zeroRow)
710 {
711 globalToLocal(row) = numZeroRows;
712 numZeroRows++;
713 }
714 }
715
716 // if more than 2 coefs are zero, continue
717
718 if (numZeroRows < 2)
719 {
720 return;
721 }
722
723 // evaluate basis functions with zero coefs on grevielle points
724
725 gsMatrix<> anchors;
726 bspline.basis().anchors_into(anchors);
727
728 gsMatrix<> params(2, coefs.rows());
729 gsVector<> localToGlobal(coefs.rows());
730 index_t counter = 0;
731
732 gsMatrix<> support = bspline.basis().support();
733
734 for (index_t fun = 0; fun != globalToLocal.rows(); fun++)
735 {
736 if (globalToLocal(fun) != -1)
737 {
738 params.col(counter) = anchors.col(fun);
739 localToGlobal(counter) = fun;
740 for (index_t row = 0; row != support.rows(); row++)
741 {
742 if (params(row, counter) < support(row, 0))
743 {
744 params(row, counter) = support(row, 0);
745 }
746
747 if (support(row, 1) < params(row, counter))
748 {
749 params(row, counter) = support(row, 1);
750 }
751 }
752 counter++;
753 }
754 }
755
756 params.conservativeResize(2, counter);
757 localToGlobal.conservativeResize(counter);
758
759 gsMatrix<> points;
760 surface.eval_into(params, points);
761
762 // do least square fit
763
764 gsSparseMatrix<> A(localToGlobal.rows(), localToGlobal.rows());
765 A.setZero();
766 gsMatrix<> B(localToGlobal.rows(), surface.geoDim());
767 B.setZero();
768
769 gsMatrix<> value;
770 gsMatrix<index_t> actives;
771
772 for (index_t k = 0; k != params.cols(); k++)
773 {
774 const gsMatrix<>& curr_param = params.col(k);
775
776 bspline.basis().eval_into(curr_param, value);
777 bspline.basis().active_into(curr_param, actives);
778 const index_t numActive = actives.rows();
779
780 for (index_t i = 0; i != numActive; i++)
781 {
782 const int I = globalToLocal(actives(i, 0));
783 if (I != -1)
784 {
785 B.row(I) += value(i, 0) * points.col(k);
786
787 for (index_t j = 0; j != numActive; j++)
788 {
789 const int J = globalToLocal(actives(j, 0));
790 if (J != -1)
791 {
792 A(I, J) += value(i, 0) * value(j, 0);
793 }
794 }
795 }
796 }
797 }
798
799 A.makeCompressed();
800 gsSparseSolver<real_t>::BiCGSTABILUT solver(A);
801 if (solver.preconditioner().info() != gsEigen::Success)
802 {
803 gsWarn<< "The preconditioner failed. Aborting.\n";
804 return;
805 }
806
807 gsMatrix<> x = solver.solve(B);
808
809 for (index_t row = 0; row != x.rows(); row++)
810 {
811 coefs.row(localToGlobal(row)) = x.row(row);
812 }
813}
814
815template <class T>
816bool exportCheck(const gsTHBSpline<2, T>& surface)
817{
818 typedef typename gsKnotVector<T>::mult_t mult_t;
819
820 for (index_t dim = 0; dim != 2; dim++)
821 {
822 typedef std::vector< gsTensorBSplineBasis< 2, real_t>* > Bases;
823 const Bases& bases = surface.basis().getBases();
824 const int deg = (bases[0])->degree(dim);
825 std::vector<mult_t> mult = (bases[0])->knots(dim).multiplicities();
826
827 if (!validMultiplicities(mult, deg))
828 {
829 return false;
830 }
831 }
832
833 return true;
834}
835
836template <class T>
837bool exportTHBsurface( const gsTHBSpline<2, T>& surface,
838 std::vector<gsTrimData<T> >& trimData,
839 PK_ASSEMBLY_t& assembly )
840{
841 const gsTHBSplineBasis<2>& basis= surface.basis();
842
843 PK_ERROR_t err = PK_SESSION_set_check_continuity(PK_LOGICAL_false);
844 PARASOLID_ERROR(PK_SESSION_set_check_continuity, err);
845
846 err = PK_ASSEMBLY_create_empty(&assembly);
847 PARASOLID_ERROR(PK_ASSEMBLY_create_empty, err);
848
849 for (size_t box = 0;box<trimData.size();++box)
850 {
851 unsigned level = trimData[box].m_level;
852 std::vector<index_t> AABBBox = trimData[box].m_AABBBox;
853 std::vector<std::vector<std::vector<T> > > polylines = trimData[box].m_Polylines;
854
855 gsTensorBSpline<2, T> bspline =
856 basis.getBSplinePatch(AABBBox, level, surface.coefs());
857
858 //std::cout << "bspline size: " << bspline.basis().size() << std::endl;
859 //std::cout << "bspline coef size: " << bspline.coefs().rows() << std::endl;
860 makeValidGeometry(surface, bspline);
861
862 //gsWriteParaview(bspline,"paraviewtest" + util::to_string(box),1000,true,true);
863
864 PK_BSURF_t bsurf;
865 createPK_BSURF<T>(bspline, bsurf); // swap
866
867 std::vector<PK_CURVE_t> curves;
868 std::vector<PK_INTERVAL_t> intervals;
869 std::vector<int> trim_loop;
870 std::vector<int> trim_set;
871
872 for (unsigned loop = 0; loop != polylines.size(); loop++)
873 {
874 for (unsigned seg = 0; seg != polylines[loop].size(); seg++)
875 {
876 real_t y1 = polylines[loop][seg][0];
877 real_t x1 = polylines[loop][seg][1];
878 real_t y2 = polylines[loop][seg][2];
879 real_t x2 = polylines[loop][seg][3];
880
881 PK_CURVE_t line;
882 PK_INTERVAL_t intervalDummy;
883 PK_INTERVAL_t interval;
884
885 PK_SURF_make_curve_isoparam_o_t options;
886 PK_SURF_make_curve_isoparam_o_m(options);
887
888 if (x1 == x2)
889 {
890 err = PK_SURF_make_curve_isoparam(bsurf, x1, PK_PARAM_direction_v_c,
891 &options, &line, &intervalDummy);
892 PARASOLID_ERROR(PK_SURF_make_curve_isoparam, err);
893
894 getInterval(true, y1, y2, x1, bspline, line, interval);
895 }
896 else
897 {
898 err = PK_SURF_make_curve_isoparam(bsurf, y1, PK_PARAM_direction_u_c,
899 &options, &line, &intervalDummy);
900 PARASOLID_ERROR(PK_SURF_make_curve_isoparam, err);
901
902 getInterval(false, x1, x2, y1, bspline, line, interval);
903 }
904
905 curves.push_back(line);
906 intervals.push_back(interval);
907 trim_loop.push_back(loop);
908 trim_set.push_back(0);
909
910 // ------------------------------------------------------------
911 // very usefull for debugging / logging
912// PK_PART_transmit_o_t transmit_options;
913// PK_PART_transmit_o_m(transmit_options);
914// transmit_options.transmit_format = PK_transmit_format_text_c;
915
916// std::string name = "level_" + internal::to_string(level) +
917// "_box_" + internal::to_string(box) +
918// "_curve_" + internal::to_string(curve) +
919// "_seg_" + internal::to_string(seg);
920
921// PK_GEOM_copy_o_t copyOptions;
922// PK_GEOM_copy_o_m(copyOptions);
923
924// PK_GEOM_copy_r_t result;
925// err = PK_GEOM_copy(1, &line, &copyOptions, &result);
926// PARASOLID_ERROR(PK_GEOM_copy, err);
927
928// PK_CURVE_t copyLine = *(result.copied_geoms);
929
930// PK_BODY_t lineBody;
931// PK_CURVE_make_wire_body_o_t option;
932// PK_CURVE_make_wire_body_o_m(option);
933// int line_new_edges;
934// PK_EDGE_t** edges = NULL;
935// int** edge_index = NULL;
936
937// err = PK_CURVE_make_wire_body_2(1, &copyLine, &interval, &option,
938// &lineBody, &line_new_edges, edges, edge_index);
939// PARASOLID_ERROR(PK_CURVE_make_wire_body_2, err);
940
941
942
943// PK_ERROR_code_t err = PK_PART_transmit(1, &lineBody, name.c_str(), &transmit_options);
944// PARASOLID_ERROR(PK_PART_transmit, err);
945 // ------------------------------------------------------------
946
947 }
948 }
949 // ----------------------------------------------------------------------
950 // very usefull for debugging / logging
951
952 // copy
953// PK_GEOM_copy_o_t copyOptions;
954// PK_GEOM_copy_o_m(copyOptions);
955// PK_GEOM_copy_r_t result;
956// err = PK_GEOM_copy(1, &bsurf, &copyOptions, &result);
957// PARASOLID_ERROR(PK_GEOM_copy, err);
958// PK_BSURF_t copyBsurf = *(result.copied_geoms);
959
960
961// PK_ASSEMBLY_t part;
962// PK_ASSEMBLY_create_empty(&part);
963// err = PK_PART_add_geoms(part, 1, &copyBsurf);
964// PARASOLID_ERROR(PK_PART_add_geoms, err);
965
966// PK_PART_transmit_o_t transmit_options;
967// PK_PART_transmit_o_m(transmit_options);
968// transmit_options.transmit_format = PK_transmit_format_text_c;
969
970
971// std::string name = "level_" + internal::to_string(level) +
972// "_box_" + internal::to_string(box);
973
974// PK_ERROR_code_t err = PK_PART_transmit(1, &part, name.c_str(), &transmit_options);
975// PARASOLID_ERROR(PK_PART_transmit, err);
976
977 // ----------------------------------------------------------------------
978
979 // Extra transposition (2019-11-14) to compensate for the transposition in createPK_BSURF.
980 // PK_BSURF_reparameterise_o_t rep_options;
981 // PK_BSURF_reparameterise_o_m(rep_options);
982 // rep_options.transpose = PK_LOGICAL_true;
983
984 // err = PK_BSURF_reparameterise(bsurf,&rep_options);
985 // PARASOLID_ERROR(PK_BSURF_reparameterise, err);
986 // End of the transposition.
987
988
989 PK_SURF_trim_data_t trim_data;
990 trim_data.n_spcurves = static_cast<int>(curves.size());
991 trim_data.spcurves = curves.data();
992 trim_data.intervals = intervals.data();
993 trim_data.trim_loop = trim_loop.data();
994 trim_data.trim_set = trim_set.data();
995
996 PK_SURF_make_sheet_trimmed_o_t options;
997 PK_SURF_make_sheet_trimmed_o_m(options);
998 options.check_loops = PK_LOGICAL_true;
999
1000 PK_BODY_t trimSurface;
1001 PK_check_state_t state;
1002 err = PK_SURF_make_sheet_trimmed(bsurf, trim_data, 1e-8, &options,
1003 &trimSurface, &state);
1004 PARASOLID_ERROR(PK_SURF_make_sheet_trimmed, err);
1005 // set the name of the body to identify it
1006 std::string name = "BoxID_" + internal::to_string(box);
1007 PK_ATTDEF_t attdef;
1008 PK_ATTRIB_t attribut;
1009 PK_ATTDEF_find("SDL/TYSA_NAME",&attdef);
1010 if(attdef == PK_ENTITY_null)
1011 std::cout<<"entity null for tysa name" <<std::endl;
1012 PK_ERROR_code_t err2 = PK_ATTRIB_create_empty(trimSurface,attdef,&attribut);
1013 if(err2 == PK_ERROR_existing_attrib || err2 == PK_ERROR_wrong_entity)
1014 std::cout<<"attribute already exists" <<std::endl;
1015 PK_ATTRIB_set_string(attribut,0,name.c_str());
1016 if (state != PK_BODY_state_ok_c)
1017 {
1018 std::cout << "Something went wrong. state(" << state << ")" << std::endl;
1019 }
1020
1021 PK_INSTANCE_sf_t sform;
1022 sform.assembly = assembly;
1023 sform.transf = PK_ENTITY_null;
1024 sform.part = trimSurface;
1025
1026 PK_INSTANCE_t instance;
1027 err = PK_INSTANCE_create(&sform, &instance);
1028 PARASOLID_ERROR(PK_INSTANCE_create, err);
1029
1030 }
1031
1032 err = PK_SESSION_set_check_continuity(PK_LOGICAL_true);
1033 PARASOLID_ERROR(PK_SESSION_set_check_continuity, err);
1034
1035 return true;
1036}
1037
1038template <class T>
1039bool getParBoxAsIndexBoxInLevel(const gsTHBSplineBasis<2, T>& basis,unsigned lvl,const std::vector<real_t>& par_box,
1040 std::vector<index_t>& index_box)
1041{
1042 T lowU=par_box[0];
1043 T lowV=par_box[1];
1044 T upU=par_box[2];
1045 T upV=par_box[3];
1046 unsigned lowerIndexU,upperIndexU,lowerIndexV,upperIndexV;
1047 const typename gsBSplineTraits<2,T>::Basis & tBasis = *(basis.getBases()[lvl]);
1048 const gsKnotVector<T> & tKvU = tBasis.component(0).knots();
1049 if( !( *(tKvU.begin())<=lowU && lowU<=upU &&
1050 upU<=tKvU.get().end()[-tKvU.degree()-1]) )
1051 {
1052 if( !(lowU<=upU) )
1053 {
1054 std::cout << "Error in the box." << std::endl;
1055 return false;
1056 }
1057 std::cout<<"Box in u is not in Range, it will be cut at domain boundary."<<std::endl;
1058 if( !( *(tKvU.begin())<=lowU) )
1059 lowU=*(tKvU.begin());
1060 if( !(upU<=tKvU.get().end()[-tKvU.degree()-1]) )
1061 upU=tKvU.get().end()[-tKvU.degree()-1];
1062 }
1063 lowerIndexU=(tKvU.uFind(lowU)-tKvU.uFind(0));
1064 upperIndexU=(tKvU.uFind(upU)-tKvU.uFind(0))+1;
1065
1066 const gsKnotVector<T> & tKvV = tBasis.component(1).knots();
1067 if( !( *(tKvV.begin())<=lowV && lowV<=upV &&
1068 upV<=tKvV.get().end()[-tKvV.degree()-1]) )
1069 {
1070 if( !(lowV<=upV) )
1071 {
1072 std::cout << "Error in the box." << std::endl;
1073 return false;
1074 }
1075 std::cout<<"Box in v is not in Range, it will be cut at domain boundary."<<std::endl;
1076 if( !( *(tKvV.begin())<=lowV) )
1077 lowV=*(tKvV.begin());
1078 if( !(upV<=tKvV.get().end()[-tKvV.degree()-1]) )
1079 upV=tKvV.get().end()[-tKvV.degree()-1];
1080 }
1081 lowerIndexV=(tKvV.uFind(lowV)-tKvV.uFind(0));
1082 upperIndexV=(tKvV.uFind(upV)-tKvV.uFind(0))+1;
1083 index_box.clear();
1084 index_box.push_back(lvl);
1085 index_box.push_back(lowerIndexU);
1086 index_box.push_back(lowerIndexV);
1087 index_box.push_back(upperIndexU);
1088 index_box.push_back(upperIndexV);
1089 return true;
1090}
1091
1092template <class T>
1093bool parBoxesIntersect(const std::vector<T>& par_boxes)
1094{
1095 T lowerUi,lowerVi,upperUi,upperVi,lowerUj,lowerVj,upperUj,upperVj;
1096 unsigned d=2;
1097 unsigned boxSize=2*d;
1098 for(unsigned i = 0;i<par_boxes.size();i+=boxSize)
1099 {
1100 lowerUi=par_boxes[i];
1101 lowerVi=par_boxes[i+1];
1102 upperUi=par_boxes[i+2];
1103 upperVi=par_boxes[i+3];
1104 for(unsigned j = 0;j<i;j+=boxSize)
1105 {
1106 lowerUj=par_boxes[j];
1107 lowerVj=par_boxes[j+1];
1108 upperUj=par_boxes[j+2];
1109 upperVj=par_boxes[j+3];
1110
1111 if(upperUi>lowerUj&&lowerUi<upperUj&&
1112 upperVi>lowerVj&&lowerVi<upperVj)
1113 return true;
1114 }
1115 }
1116 return false;
1117}
1118
1119template <class T>
1120bool getTrimDataFromBoundaryTrimCurves(gsTHBSplineBasis<2>::AxisAlignedBoundingBox boundaryAABB,
1121gsTHBSplineBasis<2>::TrimmingCurves trimCurves,std::vector<gsTrimData<T> >& trimdata)
1122{
1123 trimdata.clear();
1124 for(unsigned level = 0;level<trimCurves.size();++level)
1125 {
1126 for(unsigned component=0;component<trimCurves[level].size();++component)
1127 {
1128 gsTrimData<T> td(level,boundaryAABB[level][component],trimCurves[level][component]);
1129 trimdata.push_back(td);
1130 }
1131 }
1132 return true;
1133}
1134
1135template <class T>
1136bool getTrimCurvesAndBoundingBoxes(const gsTHBSpline<2, T>& surface,
1137 const std::vector<T>& par_boxes,
1138 std::vector<gsTrimData<T> >& trimdata)
1139{
1140 if(parBoxesIntersect(par_boxes))
1141 {
1142 std::cout<<"Given boxes have intersections."<<std::endl;
1143 return false;
1144 }
1145
1146 const gsTHBSplineBasis<2, T>* basis = static_cast< const gsTHBSplineBasis<2,T>* > (&surface.basis());
1147 unsigned maxLevel = basis->tree().getMaxInsLevel();
1148 unsigned lvl;
1149 std::vector<std::vector<index_t> >aabbBoxesForCheck;
1150 std::vector<T> par_box;
1151 std::vector<index_t> index_box;
1152 unsigned d=2;
1153 unsigned boxSize=2*d;
1154 bool success;
1155 trimdata.clear();
1156 for(unsigned i = 0;i<par_boxes.size();i=i+boxSize)
1157 {
1158 par_box.clear();
1159 par_box.push_back(par_boxes[i]);
1160 par_box.push_back(par_boxes[i+1]);
1161 par_box.push_back(par_boxes[i+2]);
1162 par_box.push_back(par_boxes[i+3]);
1163 success=getParBoxAsIndexBoxInLevel(*basis,maxLevel,par_box,index_box);
1164
1165 if(!success)
1166 return false;
1167
1168 gsVector<index_t,2>lower;
1169 lower << index_box[1],index_box[2];
1170 gsVector<index_t,2>upper;
1171 upper << index_box[3],index_box[4];
1172 lvl=basis->tree().query4(lower, upper, index_box[0]);
1173
1174 success=getParBoxAsIndexBoxInLevel(*basis,lvl,par_box,index_box);
1175 if(!success)
1176 return false;
1177
1178 //std::cout << "par-box: "<< par_box[0] << " "<< par_box[1] << " " << par_box[2] << " " << par_box[3] << std::endl;
1179 //std::cout << "el-box: "<< index_box[1] << " "<< index_box[2] << " "
1180 // << index_box[3] << " " << index_box[4] << " in level: " << lvl << std::endl;
1181
1182 std::vector<index_t> aabb_box;
1183 aabb_box.push_back(index_box[1]<<(maxLevel-lvl));
1184 aabb_box.push_back(index_box[2]<<(maxLevel-lvl));
1185 aabb_box.push_back(index_box[3]<<(maxLevel-lvl));
1186 aabb_box.push_back(index_box[4]<<(maxLevel-lvl));
1187
1188 std::vector<std::vector<std::vector<T> > > trimCurveComp;
1189 std::vector<std::vector<T> >trimCurveBox;
1190 std::vector<T> trimCurve;
1191 trimCurve.push_back(par_box[0]);
1192 trimCurve.push_back(par_box[1]);
1193 trimCurve.push_back(par_box[0]);
1194 trimCurve.push_back(par_box[3]);
1195 trimCurveBox.push_back(trimCurve);
1196 trimCurve.clear();
1197 trimCurve.push_back(par_box[0]);
1198 trimCurve.push_back(par_box[3]);
1199 trimCurve.push_back(par_box[2]);
1200 trimCurve.push_back(par_box[3]);
1201 trimCurveBox.push_back(trimCurve);
1202 trimCurve.clear();
1203 trimCurve.push_back(par_box[2]);
1204 trimCurve.push_back(par_box[1]);
1205 trimCurve.push_back(par_box[2]);
1206 trimCurve.push_back(par_box[3]);
1207 trimCurveBox.push_back(trimCurve);
1208 trimCurve.clear();
1209 trimCurve.push_back(par_box[0]);
1210 trimCurve.push_back(par_box[1]);
1211 trimCurve.push_back(par_box[2]);
1212 trimCurve.push_back(par_box[1]);
1213 trimCurveBox.push_back(trimCurve);
1214 trimCurveComp.push_back(trimCurveBox);
1215
1216 gsTrimData<T> td(lvl,aabb_box,trimCurveComp);
1217 trimdata.push_back(td);
1218
1219 aabbBoxesForCheck.push_back(aabb_box);
1220 }
1221 return true;
1222}
1223
1224
1225
1226}//extensions
1227
1228}//gismo
1229
1230#undef PARASOLID_ERROR
gismo::gsBSpline
A B-spline function of one argument, with arbitrary target dimension.
Definition gsBSpline.h:51
gismo::gsBSpline::degree
short_t degree(short_t i=0) const
Returns the degree of the B-spline.
Definition gsBSpline.h:186
gismo::gsGeometry
Abstract base class representing a geometry map.
Definition gsGeometry.h:93
gismo::gsGeometry::coefs
gsMatrix< T > & coefs()
Definition gsGeometry.h:340
gismo::gsGeometry::setCoefs
void setCoefs(gsMatrix< T > cc)
Set the coefficient matrix of the geometry, taking ownership of the matrix.
Definition gsGeometry.h:368
gismo::gsGeometry::geoDim
short_t geoDim() const
Dimension n of the absent physical space.
Definition gsGeometry.h:292
gismo::gsGeometry::basis
virtual const gsBasis< T > & basis() const =0
Returns a const reference to the basis of the geometry.
gismo::gsGeometry::parameterRange
gsMatrix< T > parameterRange() const
Returns the range of parameters as a matrix with two columns, [lower upper].
Definition gsGeometry.hpp:198
gismo::gsHTensorBasis::size
index_t size() const
The number of basis functions in this basis.
Definition gsHTensorBasis.hpp:298
gismo::gsKnotVector
Class for representing a knot vector.
Definition gsKnotVector.h:80
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMesh
Class Representing a triangle mesh with 3D vertices.
Definition gsMesh.h:32
gismo::gsTHBSplineBasis
Truncated hierarchical B-spline basis.
Definition gsTHBSplineBasis.h:36
gismo::gsTHBSpline
A truncated hierarchical B-Spline function, in d dimensions.
Definition gsTHBSpline.h:38
gismo::gsTensorBSpline
A tensor product of d B-spline functions, with arbitrary target dimension.
Definition gsTensorBSpline.h:45
gsBSplineBasis.h
Provides declaration of BSplineBasis class.
gsBSpline.h
Represents a B-spline curve/function with one parameter.
index_t
#define index_t
Definition gsConfig.h:32
gsWarn
#define gsWarn
Definition gsDebug.h:50
gsInfo
#define gsInfo
Definition gsDebug.h:43
gsFrustrum.h
Defines the Parasolud frustrim.
gsKnotVector.h
Knot vector for B-splines.
gsMesh.h
Provides declaration of the Mesh class.
gsMultiPatch.h
Provides declaration of the MultiPatch class.
gsPKSession.h
Manages starting and stopping Parasolid session.
gsReadParasolid.h
Provides declaration of gsReadParasolid functions.
gsSurface.h
Provides declaration of Surface abstract interface.
gsTHBSplineBasis.h
Provides declaration of THBSplineBasis class.
gsTHBSpline.h
Provides declaration of THBSplineBasis class.
gsTensorBSpline.h
Represents a tensor-product B-spline patch.
gismo::extensions::gsWriteParasolid
bool gsWriteParasolid(const gsGeometry< T > &ggeo, std::string const &filename)
Definition gsWriteParasolid.hpp:178
gismo::extensions::gsWritePK_SHEET
bool gsWritePK_SHEET(const gsTensorBSpline< 2, T > &tp, const std::string &filename)
Converts tp into a PK_SHEET and writes it to filename.xmt_txt.
Definition gsWriteParasolid.hpp:284
gismo::extensions::createPK_GEOM
bool createPK_GEOM(const gsGeometry< T > &ggeo, PK_GEOM_t &pgeo)
Definition gsWriteParasolid.hpp:326
gsWriteParaview.h
Provides declaration of functions writing Paraview files.
gismo::internal::to_string
std::string to_string(const unsigned &i)
Helper to convert small unsigned to string.
Definition gsXml.cpp:74
gismo
The G+Smo namespace, containing all definitions for the library.