Solvespace

1
//-----------------------------------------------------------------------------
2
// Anything involving NURBS shells (i.e., shells); except
3
// for the real math, which is in ratpoly.cpp.
4
//
5
// Copyright 2008-2013 Jonathan Westhues.
6
//-----------------------------------------------------------------------------
7
#include "../solvespace.h"
8

9
typedef struct {
10
    hSCurve     hc;
11
    hSSurface   hs;
12
} TrimLine;
13

14

15
void SShell::MakeFromExtrusionOf(SBezierLoopSet *sbls, Vector t0, Vector t1, RgbaColor color)
16
{
17
    // Make the extrusion direction consistent with respect to the normal
18
    // of the sketch we're extruding.
19
    if((t0.Minus(t1)).Dot(sbls->normal) < 0) {
20
        swap(t0, t1);
21
    }
22

23
    // Define a coordinate system to contain the original sketch, and get
24
    // a bounding box in that csys
25
    Vector n = sbls->normal.ScaledBy(-1);
26
    Vector u = n.Normal(0), v = n.Normal(1);
27
    Vector orig = sbls->point;
28
    double umax = VERY_NEGATIVE, umin = VERY_POSITIVE;
29
    sbls->GetBoundingProjd(u, orig, &umin, &umax);
30
    double vmax = VERY_NEGATIVE, vmin = VERY_POSITIVE;
31
    sbls->GetBoundingProjd(v, orig, &vmin, &vmax);
32
    // and now fix things up so that all u and v lie between 0 and 1
33
    orig = orig.Plus(u.ScaledBy(umin));
34
    orig = orig.Plus(v.ScaledBy(vmin));
35
    u = u.ScaledBy(umax - umin);
36
    v = v.ScaledBy(vmax - vmin);
37

38
    // So we can now generate the top and bottom surfaces of the extrusion,
39
    // planes within a translated (and maybe mirrored) version of that csys.
40
    SSurface s0, s1;
41
    s0 = SSurface::FromPlane(orig.Plus(t0), u, v);
42
    s0.color = color;
43
    s1 = SSurface::FromPlane(orig.Plus(t1).Plus(u), u.ScaledBy(-1), v);
44
    s1.color = color;
45
    hSSurface hs0 = surface.AddAndAssignId(&s0),
46
              hs1 = surface.AddAndAssignId(&s1);
47

48
    // Now go through the input curves. For each one, generate its surface
49
    // of extrusion, its two translated trim curves, and one trim line. We
50
    // go through by loops so that we can assign the lines correctly.
51
    SBezierLoop *sbl;
52
    for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
53
        SBezier *sb;
54
        List<TrimLine> trimLines = {};
55

56
        for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
57
            // Generate the surface of extrusion of this curve, and add
58
            // it to the list
59
            SSurface ss = SSurface::FromExtrusionOf(sb, t0, t1);
60
            ss.color = color;
61
            hSSurface hsext = surface.AddAndAssignId(&ss);
62

63
            // Translate the curve by t0 and t1 to produce two trim curves
64
            SCurve sc = {};
65
            sc.isExact = true;
66
            sc.exact = sb->TransformedBy(t0, Quaternion::IDENTITY, 1.0);
67
            (sc.exact).MakePwlInto(&(sc.pts));
68
            sc.surfA = hs0;
69
            sc.surfB = hsext;
70
            hSCurve hc0 = curve.AddAndAssignId(&sc);
71

72
            sc = {};
73
            sc.isExact = true;
74
            sc.exact = sb->TransformedBy(t1, Quaternion::IDENTITY, 1.0);
75
            (sc.exact).MakePwlInto(&(sc.pts));
76
            sc.surfA = hs1;
77
            sc.surfB = hsext;
78
            hSCurve hc1 = curve.AddAndAssignId(&sc);
79

80
            STrimBy stb0, stb1;
81
            // The translated curves trim the flat top and bottom surfaces.
82
            stb0 = STrimBy::EntireCurve(this, hc0, /*backwards=*/false);
83
            stb1 = STrimBy::EntireCurve(this, hc1, /*backwards=*/true);
84
            (surface.FindById(hs0))->trim.Add(&stb0);
85
            (surface.FindById(hs1))->trim.Add(&stb1);
86

87
            // The translated curves also trim the surface of extrusion.
88
            stb0 = STrimBy::EntireCurve(this, hc0, /*backwards=*/true);
89
            stb1 = STrimBy::EntireCurve(this, hc1, /*backwards=*/false);
90
            (surface.FindById(hsext))->trim.Add(&stb0);
91
            (surface.FindById(hsext))->trim.Add(&stb1);
92

93
            // And form the trim line
94
            Vector pt = sb->Finish();
95
            sc = {};
96
            sc.isExact = true;
97
            sc.exact = SBezier::From(pt.Plus(t0), pt.Plus(t1));
98
            (sc.exact).MakePwlInto(&(sc.pts));
99
            hSCurve hl = curve.AddAndAssignId(&sc);
100
            // save this for later
101
            TrimLine tl;
102
            tl.hc = hl;
103
            tl.hs = hsext;
104
            trimLines.Add(&tl);
105
        }
106

107
        int i;
108
        for(i = 0; i < trimLines.n; i++) {
109
            TrimLine *tl = &(trimLines[i]);
110
            SSurface *ss = surface.FindById(tl->hs);
111

112
            TrimLine *tlp = &(trimLines[WRAP(i-1, trimLines.n)]);
113

114
            STrimBy stb;
115
            stb = STrimBy::EntireCurve(this, tl->hc, /*backwards=*/true);
116
            ss->trim.Add(&stb);
117
            stb = STrimBy::EntireCurve(this, tlp->hc, /*backwards=*/false);
118
            ss->trim.Add(&stb);
119

120
            (curve.FindById(tl->hc))->surfA = ss->h;
121
            (curve.FindById(tlp->hc))->surfB = ss->h;
122
        }
123
        trimLines.Clear();
124
    }
125
}
126

127
bool SShell::CheckNormalAxisRelationship(SBezierLoopSet *sbls, Vector pt, Vector axis, double da, double dx)
128
// Check that the direction of revolution/extrusion ends up parallel to the normal of
129
// the sketch, on the side of the axis where the sketch is.
130
{
131
    SBezierLoop *sbl;
132
    Vector pto;
133
    double md = VERY_NEGATIVE;
134
    for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
135
        SBezier *sb;
136
        for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
137
            // Choose the point farthest from the axis; we'll get garbage
138
            // if we choose a point that lies on the axis, for example.
139
            // (And our surface will be self-intersecting if the sketch
140
            // spans the axis, so don't worry about that.)
141
            for(int i = 0; i <= sb->deg; i++) {
142
                Vector p = sb->ctrl[i];
143
                double d = p.DistanceToLine(pt, axis);
144
                if(d > md) {
145
                    md  = d;
146
                    pto = p;
147
                }
148
            }
149
        }
150
    }
151
    Vector ptc = pto.ClosestPointOnLine(pt, axis),
152
           up = axis.Cross(pto.Minus(ptc)).ScaledBy(da),
153
           vp = up.Plus(axis.ScaledBy(dx));
154
   
155
    return (vp.Dot(sbls->normal) > 0);
156
}
157

158
// sketch must not contain the axis of revolution as a non-construction line for helix
159
void SShell::MakeFromHelicalRevolutionOf(SBezierLoopSet *sbls, Vector pt, Vector axis,
160
                                         RgbaColor color, Group *group, double angles,
161
                                         double anglef, double dists, double distf) {
162
    int i0 = surface.n; // number of pre-existing surfaces
163
    SBezierLoop *sbl;
164
    // for testing - hard code the axial distance, and number of sections.
165
    // distance will need to be parameters in the future.
166
    double dist  = distf - dists;
167
    int sections = (int)(fabs(anglef - angles) / (PI / 2) + 1);
168
    double wedge = (anglef - angles) / sections;
169
    int startMapping = Group::REMAP_LATHE_START, endMapping = Group::REMAP_LATHE_END;
170

171
    if(CheckNormalAxisRelationship(sbls, pt, axis, anglef-angles, distf-dists)) {
172
        swap(angles, anglef);
173
        swap(dists, distf);
174
        dist  = -dist;
175
        wedge = -wedge;
176
        swap(startMapping, endMapping);
177
    }
178

179
    // Define a coordinate system to contain the original sketch, and get
180
    // a bounding box in that csys
181
    Vector n = sbls->normal.ScaledBy(-1);
182
    Vector u = n.Normal(0), v = n.Normal(1);
183
    Vector orig = sbls->point;
184
    double umax = VERY_NEGATIVE, umin = VERY_POSITIVE;
185
    sbls->GetBoundingProjd(u, orig, &umin, &umax);
186
    double vmax = VERY_NEGATIVE, vmin = VERY_POSITIVE;
187
    sbls->GetBoundingProjd(v, orig, &vmin, &vmax);
188
    // and now fix things up so that all u and v lie between 0 and 1
189
    orig = orig.Plus(u.ScaledBy(umin));
190
    orig = orig.Plus(v.ScaledBy(vmin));
191
    u    = u.ScaledBy(umax - umin);
192
    v    = v.ScaledBy(vmax - vmin);
193

194
    // So we can now generate the end caps of the extrusion within
195
    // a translated and rotated (and maybe mirrored) version of that csys.
196
    SSurface s0, s1;
197
    s0 = SSurface::FromPlane(orig.RotatedAbout(pt, axis, angles).Plus(axis.ScaledBy(dists)),
198
                             u.RotatedAbout(axis, angles), v.RotatedAbout(axis, angles));
199
    s0.color = color;
200

201
    hEntity face0 = group->Remap(Entity::NO_ENTITY, startMapping);
202
    s0.face = face0.v;
203

204
    s1 = SSurface::FromPlane(
205
        orig.Plus(u).RotatedAbout(pt, axis, anglef).Plus(axis.ScaledBy(distf)),
206
        u.ScaledBy(-1).RotatedAbout(axis, anglef), v.RotatedAbout(axis, anglef));
207
    s1.color = color;
208

209
    hEntity face1 = group->Remap(Entity::NO_ENTITY, endMapping);
210
    s1.face = face1.v;
211

212
    hSSurface hs0 = surface.AddAndAssignId(&s0);
213
    hSSurface hs1 = surface.AddAndAssignId(&s1);
214

215
    // Now we actually build and trim the swept surfaces. One loop at a time.
216
    for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
217
        int i, j;
218
        SBezier *sb;
219
        List<std::vector<hSSurface>> hsl = {};
220

221
        // This is where all the NURBS are created and Remapped to the generating curve
222
        for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
223
            std::vector<hSSurface> revs(sections);
224
            for(j = 0; j < sections; j++) {
225
                if((dist == 0) && sb->deg == 1 &&
226
                   (sb->ctrl[0]).DistanceToLine(pt, axis) < LENGTH_EPS &&
227
                   (sb->ctrl[1]).DistanceToLine(pt, axis) < LENGTH_EPS) {
228
                    // This is a line on the axis of revolution; it does
229
                    // not contribute a surface.
230
                    revs[j].v = 0;
231
                } else {
232
                    SSurface ss = SSurface::FromRevolutionOf(
233
                        sb, pt, axis, angles + (wedge)*j, angles + (wedge) * (j + 1),
234
                        dists + j * dist / sections, dists + (j + 1) * dist / sections);
235
                    ss.color = color;
236
                    if(sb->entity != 0) {
237
                        hEntity he;
238
                        he.v          = sb->entity;
239
                        hEntity hface = group->Remap(he, Group::REMAP_LINE_TO_FACE);
240
                        if(SK.entity.FindByIdNoOops(hface) != NULL) {
241
                            ss.face = hface.v;
242
                        }
243
                    }
244
                    revs[j] = surface.AddAndAssignId(&ss);
245
                }
246
            }
247
            hsl.Add(&revs);
248
        }
249
        // Still the same loop. Need to create trim curves
250
        for(i = 0; i < sbl->l.n; i++) {
251
            std::vector<hSSurface> revs = hsl[i], revsp = hsl[WRAP(i - 1, sbl->l.n)];
252

253
            sb = &(sbl->l[i]);
254

255
            // we will need the grid t-values for this entire row of surfaces
256
            List<double> t_values;
257
            t_values = {};
258
            if (revs[0].v) { 
259
                double ps = 0.0;
260
                t_values.Add(&ps);
261
                (surface.FindById(revs[0]))->MakeTriangulationGridInto(
262
                        &t_values, 0.0, 1.0, true, 0);
263
            }
264
            // we generate one more curve than we did surfaces
265
            for(j = 0; j <= sections; j++) {
266
                SCurve sc;
267
                Quaternion qs = Quaternion::From(axis, angles + wedge * j);
268
                // we want Q*(x - p) + p = Q*x + (p - Q*p)
269
                Vector ts =
270
                    pt.Minus(qs.Rotate(pt)).Plus(axis.ScaledBy(dists + j * dist / sections));
271

272
                // If this input curve generated a surface, then trim that
273
                // surface with the rotated version of the input curve.
274
                if(revs[0].v) { // not d[j] because crash on j==sections
275
                    sc         = {};
276
                    sc.isExact = true;
277
                    sc.exact   = sb->TransformedBy(ts, qs, 1.0);
278
                    // make the PWL for the curve based on t value list
279
                    for(int x = 0; x < t_values.n; x++) {
280
                        SCurvePt scpt;
281
                        scpt.tag    = 0;
282
                        scpt.p      = sc.exact.PointAt(t_values[x]);
283
                        scpt.vertex = (x == 0) || (x == (t_values.n - 1));
284
                        sc.pts.Add(&scpt);
285
                    }
286

287
                    // the surfaces already exists so trim with this curve
288
                    if(j < sections) {
289
                        sc.surfA = revs[j];
290
                    } else {
291
                        sc.surfA = hs1; // end cap
292
                    }
293

294
                    if(j > 0) {
295
                        sc.surfB = revs[j - 1];
296
                    } else {
297
                        sc.surfB = hs0; // staring cap
298
                    }
299

300
                    hSCurve hcb = curve.AddAndAssignId(&sc);
301

302
                    STrimBy stb;
303
                    stb = STrimBy::EntireCurve(this, hcb, /*backwards=*/true);
304
                    (surface.FindById(sc.surfA))->trim.Add(&stb);
305
                    stb = STrimBy::EntireCurve(this, hcb, /*backwards=*/false);
306
                    (surface.FindById(sc.surfB))->trim.Add(&stb);
307
                } else if(j == 0) { // curve was on the rotation axis and is shared by the end caps.
308
                    sc         = {};
309
                    sc.isExact = true;
310
                    sc.exact   = sb->TransformedBy(ts, qs, 1.0);
311
                    (sc.exact).MakePwlInto(&(sc.pts));
312
                    sc.surfA    = hs1; // end cap
313
                    sc.surfB    = hs0; // staring cap
314
                    hSCurve hcb = curve.AddAndAssignId(&sc);
315

316
                    STrimBy stb;
317
                    stb = STrimBy::EntireCurve(this, hcb, /*backwards=*/true);
318
                    (surface.FindById(sc.surfA))->trim.Add(&stb);
319
                    stb = STrimBy::EntireCurve(this, hcb, /*backwards=*/false);
320
                    (surface.FindById(sc.surfB))->trim.Add(&stb);
321
                }
322

323
                // And if this input curve and the one after it both generated
324
                // surfaces, then trim both of those by the appropriate
325
                // curve based on the control points.
326
                if((j < sections) && revs[j].v && revsp[j].v) {
327
                    SSurface *ss = surface.FindById(revs[j]);
328

329
                    sc         = {};
330
                    sc.isExact = true;
331
                    sc.exact   = SBezier::From(ss->ctrl[0][0], ss->ctrl[0][1], ss->ctrl[0][2]);
332
                    sc.exact.weight[1] = ss->weight[0][1];
333
                    double max_dt = 0.5;
334
                    if (sc.exact.deg > 1) max_dt = 0.125;
335
                    (sc.exact).MakePwlInto(&(sc.pts), 0.0, max_dt);
336
                    sc.surfA = revs[j];
337
                    sc.surfB = revsp[j];
338

339
                    hSCurve hcc = curve.AddAndAssignId(&sc);
340

341
                    STrimBy stb;
342
                    stb = STrimBy::EntireCurve(this, hcc, /*backwards=*/false);
343
                    (surface.FindById(sc.surfA))->trim.Add(&stb);
344
                    stb = STrimBy::EntireCurve(this, hcc, /*backwards=*/true);
345
                    (surface.FindById(sc.surfB))->trim.Add(&stb);
346
                }
347
            }
348
            t_values.Clear();
349
        }
350

351
        hsl.Clear();
352
    }
353

354
    if(dist == 0) {
355
        MakeFirstOrderRevolvedSurfaces(pt, axis, i0);
356
    }
357
}
358

359
void SShell::MakeFromRevolutionOf(SBezierLoopSet *sbls, Vector pt, Vector axis, RgbaColor color,
360
                                  Group *group) {
361
    int i0 = surface.n; // number of pre-existing surfaces
362
    SBezierLoop *sbl;
363

364
    if(CheckNormalAxisRelationship(sbls, pt, axis, 1.0, 0.0)) {
365
        axis = axis.ScaledBy(-1);
366
    }
367

368
    // Now we actually build and trim the surfaces.
369
    for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
370
        int i, j;
371
        SBezier *sb;
372
        List<std::vector<hSSurface>> hsl = {};
373

374
        for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
375
            std::vector<hSSurface> revs(4);
376
            for(j = 0; j < 4; j++) {
377
                if(sb->deg == 1 &&
378
                    (sb->ctrl[0]).DistanceToLine(pt, axis) < LENGTH_EPS &&
379
                    (sb->ctrl[1]).DistanceToLine(pt, axis) < LENGTH_EPS)
380
                {
381
                    // This is a line on the axis of revolution; it does
382
                    // not contribute a surface.
383
                    revs[j].v = 0;
384
                } else {
385
                    SSurface ss = SSurface::FromRevolutionOf(sb, pt, axis, (PI / 2) * j,
386
                                                             (PI / 2) * (j + 1), 0.0, 0.0);
387
                    ss.color = color;
388
                    if(sb->entity != 0) {
389
                        hEntity he;
390
                        he.v = sb->entity;
391
                        hEntity hface = group->Remap(he, Group::REMAP_LINE_TO_FACE);
392
                        if(SK.entity.FindByIdNoOops(hface) != NULL) {
393
                            ss.face = hface.v;
394
                        }
395
                    }
396
                    revs[j] = surface.AddAndAssignId(&ss);
397
                }
398
            }
399
            hsl.Add(&revs);
400
        }
401

402
        for(i = 0; i < sbl->l.n; i++) {
403
            std::vector<hSSurface> revs  = hsl[i],
404
                     revsp = hsl[WRAP(i-1, sbl->l.n)];
405

406
            sb   = &(sbl->l[i]);
407

408
            for(j = 0; j < 4; j++) {
409
                SCurve sc;
410
                Quaternion qs = Quaternion::From(axis, (PI/2)*j);
411
                // we want Q*(x - p) + p = Q*x + (p - Q*p)
412
                Vector ts = pt.Minus(qs.Rotate(pt));
413

414
                // If this input curve generate a surface, then trim that
415
                // surface with the rotated version of the input curve.
416
                if(revs[j].v) {
417
                    sc = {};
418
                    sc.isExact = true;
419
                    sc.exact = sb->TransformedBy(ts, qs, 1.0);
420
                    (sc.exact).MakePwlInto(&(sc.pts));
421
                    sc.surfA = revs[j];
422
                    sc.surfB = revs[WRAP(j-1, 4)];
423

424
                    hSCurve hcb = curve.AddAndAssignId(&sc);
425

426
                    STrimBy stb;
427
                    stb = STrimBy::EntireCurve(this, hcb, /*backwards=*/true);
428
                    (surface.FindById(sc.surfA))->trim.Add(&stb);
429
                    stb = STrimBy::EntireCurve(this, hcb, /*backwards=*/false);
430
                    (surface.FindById(sc.surfB))->trim.Add(&stb);
431
                }
432

433
                // And if this input curve and the one after it both generated
434
                // surfaces, then trim both of those by the appropriate
435
                // circle.
436
                if(revs[j].v && revsp[j].v) {
437
                    SSurface *ss = surface.FindById(revs[j]);
438

439
                    sc = {};
440
                    sc.isExact = true;
441
                    sc.exact = SBezier::From(ss->ctrl[0][0],
442
                                             ss->ctrl[0][1],
443
                                             ss->ctrl[0][2]);
444
                    sc.exact.weight[1] = ss->weight[0][1];
445
                    (sc.exact).MakePwlInto(&(sc.pts));
446
                    sc.surfA = revs[j];
447
                    sc.surfB = revsp[j];
448

449
                    hSCurve hcc = curve.AddAndAssignId(&sc);
450

451
                    STrimBy stb;
452
                    stb = STrimBy::EntireCurve(this, hcc, /*backwards=*/false);
453
                    (surface.FindById(sc.surfA))->trim.Add(&stb);
454
                    stb = STrimBy::EntireCurve(this, hcc, /*backwards=*/true);
455
                    (surface.FindById(sc.surfB))->trim.Add(&stb);
456
                }
457
            }
458
        }
459

460
        hsl.Clear();
461
    }
462

463
    MakeFirstOrderRevolvedSurfaces(pt, axis, i0);
464
}
465

466
void SShell::MakeFirstOrderRevolvedSurfaces(Vector pt, Vector axis, int i0) {
467
    int i;
468

469
    for(i = i0; i < surface.n; i++) {
470
        SSurface *srf = &(surface[i]);
471

472
        // Revolution of a line; this is potentially a plane, which we can
473
        // rewrite to have degree (1, 1).
474
        if(srf->degm == 1 && srf->degn == 2) {
475
            // close start, far start, far finish
476
            Vector cs, fs, ff;
477
            double d0, d1;
478
            d0 = (srf->ctrl[0][0]).DistanceToLine(pt, axis);
479
            d1 = (srf->ctrl[1][0]).DistanceToLine(pt, axis);
480

481
            if(d0 > d1) {
482
                cs = srf->ctrl[1][0];
483
                fs = srf->ctrl[0][0];
484
                ff = srf->ctrl[0][2];
485
            } else {
486
                cs = srf->ctrl[0][0];
487
                fs = srf->ctrl[1][0];
488
                ff = srf->ctrl[1][2];
489
            }
490

491
            // origin close, origin far
492
            Vector oc = cs.ClosestPointOnLine(pt, axis),
493
                   of = fs.ClosestPointOnLine(pt, axis);
494

495
            if(oc.Equals(of)) {
496
                // This is a plane, not a (non-degenerate) cone.
497
                Vector oldn = srf->NormalAt(0.5, 0.5);
498

499
                Vector u = fs.Minus(of), v;
500

501
                v = (axis.Cross(u)).WithMagnitude(1);
502

503
                double vm = (ff.Minus(of)).Dot(v);
504
                v = v.ScaledBy(vm);
505

506
                srf->degm = 1;
507
                srf->degn = 1;
508
                srf->ctrl[0][0] = of;
509
                srf->ctrl[0][1] = of.Plus(u);
510
                srf->ctrl[1][0] = of.Plus(v);
511
                srf->ctrl[1][1] = of.Plus(u).Plus(v);
512
                srf->weight[0][0] = 1;
513
                srf->weight[0][1] = 1;
514
                srf->weight[1][0] = 1;
515
                srf->weight[1][1] = 1;
516

517
                if(oldn.Dot(srf->NormalAt(0.5, 0.5)) < 0) {
518
                    swap(srf->ctrl[0][0], srf->ctrl[1][0]);
519
                    swap(srf->ctrl[0][1], srf->ctrl[1][1]);
520
                }
521
                continue;
522
            }
523

524
            if(fabs(d0 - d1) < LENGTH_EPS) {
525
                // This is a cylinder; so transpose it so that we'll recognize
526
                // it as a surface of extrusion.
527
                SSurface sn = *srf;
528

529
                // Transposing u and v flips the normal, so reverse u to
530
                // flip it again and put it back where we started.
531
                sn.degm = 2;
532
                sn.degn = 1;
533
                int dm, dn;
534
                for(dm = 0; dm <= 1; dm++) {
535
                    for(dn = 0; dn <= 2; dn++) {
536
                        sn.ctrl  [dn][dm] = srf->ctrl  [1-dm][dn];
537
                        sn.weight[dn][dm] = srf->weight[1-dm][dn];
538
                    }
539
                }
540

541
                *srf = sn;
542
                continue;
543
            }
544
        }
545
    }
546
}
547

548
void SShell::MakeFromCopyOf(SShell *a) {
549
    ssassert(this != a, "Can't make from copy of self");
550
    MakeFromTransformationOf(a,
551
        Vector::From(0, 0, 0), Quaternion::IDENTITY, 1.0);
552
}
553

554
void SShell::MakeFromTransformationOf(SShell *a,
555
                                      Vector t, Quaternion q, double scale)
556
{
557
    booleanFailed = false;
558
    surface.ReserveMore(a->surface.n);
559
    for(SSurface &s : a->surface) {
560
        SSurface n;
561
        n = SSurface::FromTransformationOf(&s, t, q, scale, /*includingTrims=*/true);
562
        surface.Add(&n); // keeping the old ID
563
    }
564

565
    curve.ReserveMore(a->curve.n);
566
    for(SCurve &c : a->curve) {
567
        SCurve n;
568
        n = SCurve::FromTransformationOf(&c, t, q, scale);
569
        curve.Add(&n); // keeping the old ID
570
    }
571
}
572

573
void SShell::MakeEdgesInto(SEdgeList *sel) {
574
    for(SSurface &s : surface) {
575
        s.MakeEdgesInto(this, sel, SSurface::MakeAs::XYZ);
576
    }
577
}
578

579
void SShell::MakeSectionEdgesInto(Vector n, double d, SEdgeList *sel, SBezierList *sbl)
580
{
581
    for(SSurface &s : surface) {
582
        if(s.CoincidentWithPlane(n, d)) {
583
            s.MakeSectionEdgesInto(this, sel, sbl);
584
        }
585
    }
586
}
587

588
void SShell::TriangulateInto(SMesh *sm) {
589
#pragma omp parallel for
590
    for(int i=0; i<surface.n; i++) {
591
        SSurface *s = &surface[i];
592
        SMesh m;
593
        s->TriangulateInto(this, &m);
594
        #pragma omp critical
595
        sm->MakeFromCopyOf(&m);
596
        m.Clear();
597
    }
598
}
599

600
bool SShell::IsEmpty() const {
601
    return surface.IsEmpty();
602
}
603

604
void SShell::Clear() {
605
    for(SSurface &s : surface) {
606
        s.Clear();
607
    }
608
    surface.Clear();
609

610
    for(SCurve &c : curve) {
611
        c.Clear();
612
    }
613
    curve.Clear();
614
}
615