Coverage for dragonfly_sim/core/shape_design.py: 88%
367 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1"""Shape design variables as modular layers on an lsdo_geo FFD block.
3``DG_windtunnel_model`` only sees an ``FFDShapeParameterization``: it hands it
4the wall nodes in ``set_up_sim`` (``setup``) and takes the resulting wall
5displacement into the mesh warper (``deform_mesh``). Which design variables
6exist, their bounds, and how they move the control points is decided entirely
7by the layers given to the parameterization:
9* ``ControlPointMotions`` -- raw control-point motions (``cp_motions``).
10* ``SectionalVariables`` -- the ``SectionalShape`` presets (2D camber and
11 thickness per chordwise station; 3D twist, chord, ... per spanwise section).
12* ``WingShape`` -- a trapezoidal wing: quarter-chord (or any chord-line) sweep,
13 aspect ratio, root chord, span and taper ratio, plus chordwise thickness and
14 camber modes per spanwise station.
16A driver then reads::
18 shape = FFDShapeParameterization(ffd_shape, degree, corner_list,
19 layers=[WingShape(planform={...}, sections={...})])
20 model = DG_windtunnel_model(mesh, boundary_dict, shape, ...)
21 model.set_up_sim()
22 shape_parameters = shape.declare_design_variables()
23 ...
24 u_vec = model.evaluate(model.deform_mesh(), shape_parameters, alpha)
26and constrains geometry through ``shape.outputs()`` (layer quantities such as
27the wing's planform area) and ``shape.geometry`` (``WallGeometry``: enclosed
28measure, thickness at stations, planform stations).
30The coefficient pipeline is ``C0 -> one lsdo_geo SectionalParameterization ->
31additive layers``. lsdo_geo freezes axes and stretch origins from the points'
32values when the graph is built, so every layer's sectional operations are
33collected into a single evaluation on the constant baseline (and must share
34one principal direction); additive motions (cp motions, sectional modes) are
35linear and go on top.
37A new layer subclasses ``ShapeLayer`` and overrides what it needs.
38"""
39import numpy as np
40import csdl_alpha as csdl
41import lsdo_function_spaces as lfs
42import scipy.sparse as sps
43from scipy.optimize import brentq
45from dragonfly_sim.core.shape_parameterization import (WallFFD, SectionalShape, WallGeometry,
46 _inline_recording)
47from dragonfly_sim.utils.ffd_dv_utils import (ShapeDVSet, cp_dv_directions, build_cp_motion_dv,
48 build_sectional_dv, build_planform_dv)
49from dragonfly_sim.utils.meshwarping_utils import inner_bdry_function_from_corner_list
52class FFDShapeParameterization:
53 """An FFD block around the wall plus the design-variable layers that move it.
55 Parameters
56 ----------
57 shape : sequence of int
58 Control points per parametric direction.
59 degree : int or sequence of int
60 block_corner_list : optional
61 Two-surface block definition (``shape_parameterization.
62 ffd_corners_from_corner_list``); None wraps the wall's bounding box.
63 Also locates the wall when the model gets no explicit
64 ``mesh_inner_bdry_function`` (``inner_boundary_function``).
65 layers : sequence of ShapeLayer
66 Applied in order; their design variables are declared in this order.
67 margin, projection_newton_tol : see ``WallFFD``.
68 inner_bdry_rel_tol, inner_bdry_overset_margin : containment test of
69 ``inner_boundary_function``.
70 """
72 def __init__(self, shape, degree=2, block_corner_list=None, layers=(), margin=1e-4,
73 projection_newton_tol=1e-12, inner_bdry_rel_tol=1e-6,
74 inner_bdry_overset_margin=1e-4):
75 self.shape = tuple(int(n) for n in shape)
76 self.degree = degree
77 self.block_corner_list = block_corner_list
78 self.layers = list(layers)
79 if not self.layers:
80 raise ValueError("FFDShapeParameterization needs at least one layer")
81 names = [layer.name for layer in self.layers]
82 if len(set(names)) != len(names):
83 raise ValueError("layer names must be unique, got {}".format(names))
84 self.margin = margin
85 self.projection_newton_tol = projection_newton_tol
86 self.inner_bdry_rel_tol = inner_bdry_rel_tol
87 self.inner_bdry_overset_margin = inner_bdry_overset_margin
89 self.ffd = None
90 self._wall_facets = None
91 self._geometry = None
92 self.dvs = None
93 self._coefficients = None
94 self._wall_displacement = None
96 # ---- set-up (called by the model) -------------------------------------
97 def inner_boundary_function(self):
98 """Wall facet marker from the block corners, or None without them."""
99 if self.block_corner_list is None:
100 return None
101 return inner_bdry_function_from_corner_list(
102 self.block_corner_list, overset_margin=self.inner_bdry_overset_margin,
103 rel_tol=self.inner_bdry_rel_tol)
105 def setup(self, comm, wall_coords, wall_rows, wall_facets=None):
106 """Build the FFD block around the wall and set up the layers. Collective.
108 wall_coords : (M_global, dim) replicated wall node coordinates, in the
109 rank-block order of the mesh's wall node set; wall_rows : the rows
110 this rank owns; wall_facets : replicated wall facets (rows of wall
111 node ids), needed for ``geometry`` and by layers that measure the
112 baseline wall.
113 """
114 if self.ffd is not None:
115 raise RuntimeError("FFDShapeParameterization.setup called twice")
116 self.ffd = WallFFD(comm, wall_coords, wall_rows, self.shape, degree=self.degree,
117 ffd_block_corner_list=self.block_corner_list, margin=self.margin,
118 projection_newton_tol=self.projection_newton_tol)
119 self._wall_facets = wall_facets
120 for layer in self.layers:
121 layer.setup(self)
123 @property
124 def dim(self):
125 return self._require_setup().dim
127 @property
128 def geometry(self):
129 """``WallGeometry`` of the baseline wall, built on first use."""
130 self._require_setup()
131 if self._geometry is None:
132 if self._wall_facets is None:
133 raise RuntimeError("no wall facets were given to setup()")
134 self._geometry = WallGeometry(self.ffd, self._wall_facets)
135 return self._geometry
137 def _require_setup(self):
138 if self.ffd is None:
139 raise RuntimeError("the shape parameterization is not set up yet "
140 "(DG_windtunnel_model.set_up_sim calls setup)")
141 return self.ffd
143 # ---- design variables ---------------------------------------------------
144 def declare_design_variables(self):
145 """Create every layer's csdl design variables (with bounds), in layer
146 order. Returns ``parameter_vector()``. Call inside the recorder, after
147 ``set_up_sim``."""
148 self._require_setup()
149 if self.dvs is not None:
150 raise RuntimeError("design variables already declared")
151 self.dvs = ShapeDVSet()
152 for layer in self.layers:
153 layer.declare(self.dvs)
154 if not self.dvs.names:
155 raise ValueError("the layers declared no design variables")
156 return self.parameter_vector()
158 def parameter_vector(self):
159 """All shape design variables as one flat vector (declaration order)."""
160 if self.dvs is None:
161 raise RuntimeError("call declare_design_variables() first")
162 return self.dvs.flat_variable()
164 def __getitem__(self, name):
165 """A declared design variable by name."""
166 return self.dvs[name]
168 def print_design_variables(self, printer=print):
169 printer("Shape design variables ({} entries):".format(self.dvs.size))
170 for name, dv in zip(self.dvs.names, self.dvs.dvs):
171 line = " {}: shape {}, lower [{:.4g}, {:.4g}], upper [{:.4g}, {:.4g}]".format(
172 name, dv.shape, dv.lower.min(), dv.lower.max(), dv.upper.min(), dv.upper.max())
173 if dv.value.size == 1:
174 line += ", initial {:.6g}".format(float(dv.value.ravel()[0]))
175 printer(line)
176 for layer in self.layers:
177 for line in layer.describe():
178 printer(line)
180 # ---- the coefficient pipeline -------------------------------------------
181 def coefficients(self):
182 """Deformed FFD coefficients (built once, then cached)."""
183 if self._coefficients is not None:
184 return self._coefficients
185 if self.dvs is None:
186 raise RuntimeError("call declare_design_variables() first")
187 coefficients = self.ffd.baseline_variable()
189 sectional_layers = [layer for layer in self.layers if layer.principal_dim is not None]
190 principal_dims = {layer.principal_dim for layer in sectional_layers}
191 if len(principal_dims) > 1:
192 raise ValueError("layers with sectional operations must share one principal "
193 "direction (they form a single lsdo_geo evaluation on the "
194 "baseline), got {}".format(sorted(principal_dims)))
195 if sectional_layers:
196 sectional = SectionalShape(self.ffd, principal_dim=principal_dims.pop())
197 # The layers may build value-dependent expressions of their design
198 # variables for lsdo_geo; record them inline with it.
199 with _inline_recording():
200 for layer in sectional_layers:
201 layer.sectional(sectional)
202 coefficients = sectional.apply(coefficients)
204 for layer in self.layers:
205 coefficients = layer.additive(coefficients)
206 self._coefficients = coefficients
207 return coefficients
209 def wall_displacement(self):
210 """Replicated ``(M_global, dim)`` wall node displacement (cached)."""
211 if self._wall_displacement is None:
212 self._wall_displacement = self.ffd.wall_displacement(self.coefficients())
213 return self._wall_displacement
215 def outputs(self):
216 """Constraint quantities of the layers: ``{layer.name: {...}}`` for the
217 layers that have any (e.g. ``outputs()['wing']['area']``)."""
218 self.coefficients()
219 return {layer.name: layer.outputs() for layer in self.layers if layer.outputs()}
222class ShapeLayer:
223 """Base of the design-variable layers of an FFDShapeParameterization.
225 A layer with sectional (lsdo_geo) operations sets ``principal_dim`` and
226 registers them in ``sectional``; that call happens inside an inline
227 recording, so expressions of the design variables built there have values,
228 which lsdo_geo needs. ``additive`` adds motions that are linear in the
229 variables on top of the sectional result.
230 """
231 name = 'layer'
232 principal_dim = None
234 def setup(self, parameterization):
235 """Numpy precomputation once the FFD block exists (``parameterization.ffd``,
236 ``parameterization.geometry``)."""
237 self.parameterization = parameterization
239 def declare(self, dvs):
240 """Create this layer's design variables in ``dvs`` (a ShapeDVSet)."""
242 def sectional(self, sectional):
243 """Register sectional operations on ``sectional`` (a SectionalShape)."""
245 def additive(self, coefficients):
246 return coefficients
248 def outputs(self):
249 return {}
251 def describe(self):
252 """Lines for ``FFDShapeParameterization.print_design_variables``."""
253 return []
256class ControlPointMotions(ShapeLayer):
257 """Raw control-point motions: one design variable ``cp_motions`` of shape
258 ``ffd_shape + (n_opt,)``, with bounds from ``dv_spec`` (see ffd_dv_utils)."""
259 name = 'cp_motions'
261 def __init__(self, dv_spec, scaler_mode='bounds', name='cp_motions'):
262 self.dv_spec = dv_spec
263 self.coord_idxs = cp_dv_directions(dv_spec)
264 self.scaler_mode = scaler_mode
265 self.name = name
267 def declare(self, dvs):
268 ffd = self.parameterization.ffd
269 self.dv = build_cp_motion_dv(ffd.baseline_coefficients, self.dv_spec,
270 scaler_mode=self.scaler_mode)
271 self.variable = dvs.add(self.name, self.dv)
273 def additive(self, coefficients):
274 return self.parameterization.ffd.apply_cp_motions(coefficients, self.variable,
275 self.coord_idxs)
277 def describe(self):
278 axis_names = ('x', 'y', 'z')
279 return [" {}: direction {} ({})".format(self.name, int(d), axis_names[int(d)])
280 for d in self.coord_idxs]
283class SectionalVariables(ShapeLayer):
284 """``SectionalShape`` presets as design variables.
286 variables : {kind: spec}, kind one of ``SectionalShape.add``'s presets;
287 spec is the bounds (anything ``build_sectional_dv`` takes) or a dict
288 with keys ``bounds`` and optionally ``num_values`` (fewer than the
289 number of sections: a B-spline profile), ``pivot`` and
290 ``profile_degree``. Each variable is a delta from the baseline.
291 """
292 name = 'sectional'
294 def __init__(self, principal_dim, variables, scaler_mode='bounds', name='sectional'):
295 self.principal_dim = int(principal_dim)
296 if not variables:
297 raise ValueError("SectionalVariables needs at least one variable")
298 self.specs = {}
299 for kind, spec in variables.items():
300 if not isinstance(spec, dict):
301 spec = {'bounds': spec}
302 unknown = set(spec) - {'bounds', 'num_values', 'pivot', 'profile_degree'}
303 if unknown or 'bounds' not in spec:
304 raise ValueError("sectional variable {!r}: spec needs 'bounds' and may have "
305 "num_values/pivot/profile_degree, got {}".format(kind, sorted(spec)))
306 self.specs[kind] = spec
307 self.scaler_mode = scaler_mode
308 self.name = name
310 def declare(self, dvs):
311 num_sections = self.parameterization.ffd.shape[self.principal_dim]
312 self.variables = {}
313 for kind, spec in self.specs.items():
314 n = spec.get('num_values', num_sections)
315 self.variables[kind] = dvs.add(kind, build_sectional_dv(n, spec['bounds'],
316 scaler_mode=self.scaler_mode))
318 def sectional(self, sectional):
319 for kind, spec in self.specs.items():
320 sectional.add(kind, self.variables[kind], pivot=spec.get('pivot'),
321 profile_degree=spec.get('profile_degree', 2))
324def _bspline_basis(degree, n_coefficients, u):
325 space = lfs.BSplineSpace(num_parametric_dimensions=1, degree=min(degree, n_coefficients - 1),
326 coefficients_shape=(n_coefficients,))
327 basis = space.compute_basis_matrix(np.asarray(u, dtype=np.float64).reshape(-1, 1))
328 return basis.toarray() if sps.issparse(basis) else np.asarray(basis)
331def _invert_monotone_bspline(degree, values, target):
332 """Parametric coordinate u in [0, 1] where the 1D B-spline with control
333 values ``values`` equals ``target``."""
334 def f(u):
335 return float((_bspline_basis(degree, len(values), [u]) @ values)[0]) - target
336 lo, hi = f(0.), f(1.)
337 if lo * hi > 0.:
338 raise ValueError("target {} lies outside the block's range [{}, {}]"
339 .format(target, lo + target, hi + target))
340 return brentq(f, 0., 1., xtol=1e-15, rtol=4 * np.finfo(float).eps)
343def _bernstein(degree, xi):
344 """Bernstein polynomials of `degree` at xi: shape (len(xi), degree + 1)."""
345 from scipy.special import comb
346 xi = np.asarray(xi, dtype=np.float64)[:, None]
347 m = np.arange(degree + 1)[None, :]
348 return comb(degree, m) * xi ** m * (1. - xi) ** (degree - m)
351def _section_mode_profiles(kind, num_modes, xi):
352 """Chordwise profiles of the WingShape section modes at chord fractions xi:
353 shape (len(xi), num_modes). Thickness: the Bernstein polynomials of degree
354 num_modes - 1 (they sum to one: equal values scale the section uniformly).
355 Camber: 4 xi (1 - xi) times those (zero at the leading and trailing edge;
356 equal values give the parabolic camber line)."""
357 profiles = _bernstein(num_modes - 1, xi)
358 if kind == 'camber':
359 profiles = profiles * (4. * xi * (1. - xi))[:, None]
360 return profiles
363class WingShape(ShapeLayer):
364 """Planform and sectional design variables of a trapezoidal half wing.
366 The wall must be a straight-tapered half wing with its root on the y = 0
367 symmetry plane, x chordwise and z up, inside a lofted FFD block whose
368 parametric axes are (chord, span, vertical) -- the corner-list blocks of
369 the wing drivers.
371 planform : {name: (lower, upper)} for any of
372 ``sweep`` sweep angle [rad] of the chord line at
373 ``sweep_chord_fraction`` (0.25: quarter chord)
374 ``aspect_ratio`` full-wing AR = b^2 / S
375 ``root_chord``
376 ``span`` full-wing span b (twice the half wing's semi-span)
377 ``taper_ratio`` tip chord / root chord
378 Absolute values: each variable starts at the baseline wing's value and
379 its bounds must bracket it. A name left out stays at its baseline,
380 except the ``derived`` one.
381 sections : {name: spec} for ``thickness`` and/or ``camber``. spec is the
382 bounds, or a dict {'bounds': ..., 'modes': n} for n chordwise modes per
383 FFD spanwise section (default 1). The variable then has n_sections x n
384 entries, section-major (root first: all of the root's modes, then the
385 next section's); bounds as for ``build_sectional_dv`` -- a scalar, a
386 pair, or a pair of arrays of shape (n_sections,) (the same for every
387 mode of a section), (n_sections, n) or (n_sections * n,). Pin an entry
388 with (0., 0.). Deltas from the baseline, starting at zero, with xi the
389 local chord fraction and B_m the Bernstein polynomials of degree n - 1:
390 ``thickness`` relative change of the section's thickness (t/c),
391 sum_m tau_m B_m(xi) -- 0.1 in every mode is 10% thicker
392 ``camber`` change of camber over local chord, sum_m kappa_m
393 4 xi (1 - xi) B_m(xi), so the leading and trailing edge
394 stay put; kappa in every mode is a parabolic camber line
395 of maximum kappa
396 With one mode these are a uniform thickness change and a parabolic
397 camber line. The modes are least-squares fits on the block's chordwise
398 basis; they are exact when it is a Bezier (n_chord = degree + 1) of
399 degree >= n - 1 for thickness and >= n + 1 for camber (``describe``
400 reports the fit residual).
401 derived : which of ``aspect_ratio``, ``root_chord`` and ``span`` follows from
402 the other two and the taper ratio, s = AR c_r (1 + taper) / 4. Default:
403 ``span``, unless span is a variable -- then whichever of aspect_ratio
404 and root_chord is not (with only span a variable, name it).
405 sweep_chord_fraction : chord line the sweep is measured on; sections scale
406 about it, and its root point stays fixed.
407 semi_span : baseline semi-span; default the wall's largest y (which
408 includes any tip cap).
409 fit_tol : allowed deviation of the baseline leading/trailing edges from
410 straight lines, relative to the root chord.
412 The geometry: FFD section j, at baseline span fraction eta_j = y_j / s0,
413 moves to y_j' = eta_j s', its chord c_j' = c_r' (1 - (1 - taper') eta_j)
414 scales it by r_j = c_j'/c_j0 about the reference chord line in x and about
415 the chord plane in z (keeping t/c; the thickness modes act on top), and the
416 reference line is placed at x_ref,root + y_j' tan(sweep'). These are
417 lsdo_geo sectional stretches and translations; thickness and camber are
418 additive sectional modes. With the taper ratio fixed the map is affine, so
419 the deformed wall is exactly the trapezoid the variables describe; a taper
420 change is exact at the FFD sections and B-spline-blended between them.
421 """
422 name = 'wing'
423 principal_dim = 1
424 PLANFORM_KEYS = ('sweep', 'aspect_ratio', 'root_chord', 'span', 'taper_ratio')
425 SIZE_KEYS = ('aspect_ratio', 'root_chord', 'span')
426 SECTION_KEYS = ('thickness', 'camber')
428 def __init__(self, planform=None, sections=None, sweep_chord_fraction=0.25,
429 semi_span=None, derived=None, fit_tol=1e-3, scaler_mode='bounds', name='wing'):
430 self.planform_spec = dict(planform or {})
431 self.section_spec = {}
432 unknown = set(self.planform_spec) - set(self.PLANFORM_KEYS)
433 if unknown:
434 raise ValueError("unknown planform variable(s) {}; expected any of {}"
435 .format(sorted(unknown), list(self.PLANFORM_KEYS)))
436 for key, spec in dict(sections or {}).items():
437 if key not in self.SECTION_KEYS:
438 raise ValueError("unknown section variable {!r}; expected any of {}"
439 .format(key, list(self.SECTION_KEYS)))
440 if not isinstance(spec, dict):
441 spec = {'bounds': spec}
442 if 'bounds' not in spec or set(spec) - {'bounds', 'modes'}:
443 raise ValueError("section variable {!r}: spec needs 'bounds' and may have "
444 "'modes', got {}".format(key, sorted(spec)))
445 modes = int(spec.get('modes', 1))
446 if modes < 1:
447 raise ValueError("section variable {!r} needs at least one mode".format(key))
448 self.section_spec[key] = {'bounds': spec['bounds'], 'modes': modes}
449 if not self.planform_spec and not self.section_spec:
450 raise ValueError("WingShape needs at least one planform or section variable")
451 self.derived = self._derived_size(derived)
452 self.f = float(sweep_chord_fraction)
453 self.semi_span_arg = semi_span
454 self.fit_tol = fit_tol
455 self.scaler_mode = scaler_mode
456 self.name = name
458 def _derived_size(self, derived):
459 variables = [k for k in self.SIZE_KEYS if k in self.planform_spec]
460 if len(variables) == 3:
461 raise ValueError("aspect_ratio, root_chord and span are related through the taper "
462 "ratio; at most two of them can be design variables")
463 if derived is None:
464 if 'span' not in variables:
465 return 'span'
466 free = [k for k in ('aspect_ratio', 'root_chord') if k not in variables]
467 if len(free) == 2:
468 raise ValueError("with span a design variable, say which of aspect_ratio and "
469 "root_chord follows from it (derived=...); the other stays "
470 "at its baseline")
471 return free[0]
472 if derived not in self.SIZE_KEYS:
473 raise ValueError("derived must be one of {}, got {!r}".format(list(self.SIZE_KEYS),
474 derived))
475 if derived in variables:
476 raise ValueError("{!r} cannot be both a design variable and derived".format(derived))
477 return derived
479 # ---- baseline (numpy) -----------------------------------------------
480 def setup(self, parameterization):
481 super().setup(parameterization)
482 ffd = parameterization.ffd
483 if ffd.dim != 3:
484 raise ValueError("WingShape needs a 3D wall")
485 self.baseline = self._fit_baseline(parameterization.geometry)
486 self._setup_sections(ffd)
488 def _fit_baseline(self, geometry):
489 X = geometry.ffd.baseline_wall_coords
490 y_root, y_max = X[:, 1].min(), X[:, 1].max()
491 s0 = float(self.semi_span_arg) if self.semi_span_arg is not None else float(y_max)
492 if abs(y_root) > 1e-6 * s0:
493 raise ValueError("WingShape needs the wing root on the y = 0 symmetry plane; "
494 "the wall starts at y = {}".format(y_root))
495 stations = np.linspace(0.05, 0.9, 12) * s0
496 edges = [geometry.section_edges(y) for y in stations]
497 le = np.array([e[0] for e in edges])
498 te = np.array([e[1] for e in edges])
499 le_fit = np.polyfit(stations, le[:, 0], 1)
500 te_fit = np.polyfit(stations, te[:, 0], 1)
501 z_fit = np.polyfit(stations, 0.5 * (le[:, 2] + te[:, 2]), 1)
502 c_root = te_fit[1] - le_fit[1]
503 residual = max(np.abs(np.polyval(le_fit, stations) - le[:, 0]).max(),
504 np.abs(np.polyval(te_fit, stations) - te[:, 0]).max())
505 if residual > self.fit_tol * c_root:
506 raise ValueError("the wing's leading/trailing edges deviate {:.3g} from straight "
507 "lines (fit_tol {} x root chord {:.4g}); WingShape needs a "
508 "straight-tapered planform".format(residual, self.fit_tol, c_root))
509 chord_slope = te_fit[0] - le_fit[0]
510 c_tip = c_root + chord_slope * s0
511 x_ref_root = le_fit[1] + self.f * c_root
512 ref_slope = le_fit[0] + self.f * chord_slope
513 return {'sweep': float(np.arctan(ref_slope)),
514 'aspect_ratio': float(4. * s0 / (c_root + c_tip)),
515 'root_chord': float(c_root),
516 'span': 2. * s0,
517 'taper_ratio': float(c_tip / c_root),
518 'semi_span': s0,
519 'area': float((c_root + c_tip) * s0),
520 'x_ref_root': float(x_ref_root),
521 'ref_slope': float(ref_slope),
522 'x_le': le_fit, 'z_chord_plane': z_fit,
523 'fit_residual': float(residual)}
525 def _setup_sections(self, ffd):
526 C = ffd.baseline_coefficients
527 n0, n1, n2 = ffd.shape
528 deg0 = ffd.degree[0]
529 b = self.baseline
530 y = C[0, :, 0, 1]
531 tol = 1e-9 * np.abs(C).max()
532 if (np.abs(C[..., 1] - y[None, :, None]).max() > tol
533 or np.abs(C[..., 0] - C[:, :, :1, 0]).max() > tol
534 or np.abs(C[..., 2] - C[:1, :, :, 2]).max() > tol):
535 raise ValueError("WingShape needs a lofted block with parametric axes (chord, span, "
536 "vertical): spanwise sections in planes y = const, x varying "
537 "along axis 0 only and z along axis 2 only")
539 s0, cr0, lam0 = b['semi_span'], b['root_chord'], b['taper_ratio']
540 self.y = y
541 self.eta = y / s0
542 self.c0 = cr0 * (1. - (1. - lam0) * self.eta)
543 self.x_ref0 = b['x_ref_root'] + b['ref_slope'] * y
544 x_le0 = np.polyval(b['x_le'], y)
545 z_ref0 = np.polyval(b['z_chord_plane'], y)
546 # lsdo_geo normalizes a stretch by the section's control-point extent
547 self.extent_x = C[:, :, 0, 0].max(axis=0) - C[:, :, 0, 0].min(axis=0)
548 self.extent_z = C[0, :, :, 2].max(axis=1) - C[0, :, :, 2].min(axis=1)
550 # Pivots as parametric coordinates within each section: (axis 0, axis 2).
551 # lsdo_geo evaluates them on a degree-1 B-spline through the section's
552 # control points (SectionalParameterization.helpful_section_b_spline_space),
553 # not on the FFD's own basis, so invert that one.
554 self.pivots = np.zeros((n1, 2))
555 for j in range(n1):
556 self.pivots[j, 0] = _invert_monotone_bspline(1, C[:, j, 0, 0], self.x_ref0[j])
557 self.pivots[j, 1] = _invert_monotone_bspline(1, C[0, j, :, 2], z_ref0[j])
559 # Section modes: chordwise control values whose B-spline is the mode's
560 # profile in the wing chord fraction xi (least squares; exact when the
561 # chordwise basis contains it). Thickness weights also carry the control
562 # points' height above the chord plane, the distance they scale.
563 u = np.linspace(0., 1., 8 * n0 + 1)
564 basis0 = _bspline_basis(deg0, n0, u)
565 self.mode_weights = {}
566 self.mode_fit_residual = 0.
567 for key, spec in self.section_spec.items():
568 weights = np.zeros((spec['modes'],) + ffd.shape)
569 for j in range(n1):
570 xi = (basis0 @ C[:, j, 0, 0] - x_le0[j]) / self.c0[j]
571 profiles = _section_mode_profiles(key, spec['modes'], xi)
572 fit = np.linalg.lstsq(basis0, profiles, rcond=None)[0] # (n0, modes)
573 self.mode_fit_residual = max(self.mode_fit_residual,
574 float(np.abs(basis0 @ fit - profiles).max()))
575 if key == 'thickness':
576 weights[:, :, j, :] = fit.T[:, :, None] * (C[0, j, :, 2] - z_ref0[j])[None, None, :]
577 else:
578 weights[:, :, j, :] = fit.T[:, :, None]
579 self.mode_weights[key] = weights
581 # ---- design variables ---------------------------------------------------
582 def _section_bounds(self, bounds, n1, modes):
583 """Broadcast per-section (n1,) or (n1, modes) bound arrays to the
584 variable's flat (n1 * modes,) layout; scalars pass through."""
585 if not isinstance(bounds, (tuple, list)) or len(bounds) != 2:
586 return bounds
587 resolved = []
588 for bound in bounds:
589 bound = np.asarray(bound, dtype=np.float64)
590 if bound.ndim == 0:
591 resolved.append(float(bound))
592 elif bound.shape == (n1,):
593 resolved.append(np.repeat(bound, modes))
594 elif bound.size == n1 * modes:
595 resolved.append(bound.reshape(-1))
596 else:
597 raise ValueError("section bounds must be scalars or arrays of shape ({0},), "
598 "({0}, {1}) or ({2},), got {3}".format(n1, modes, n1 * modes,
599 bound.shape))
600 return tuple(resolved)
602 def declare(self, dvs):
603 self.variables = {}
604 for key in self.PLANFORM_KEYS:
605 if key in self.planform_spec:
606 dv = build_planform_dv(self.baseline[key], self.planform_spec[key],
607 scaler_mode=self.scaler_mode, name=key)
608 self.variables[key] = dvs.add(key, dv)
609 n1 = self.parameterization.ffd.shape[1]
610 for key in self.SECTION_KEYS:
611 if key in self.section_spec:
612 modes = self.section_spec[key]['modes']
613 bounds = self._section_bounds(self.section_spec[key]['bounds'], n1, modes)
614 dv = build_sectional_dv(n1 * modes, bounds, scaler_mode=self.scaler_mode)
615 self.variables[key] = dvs.add(key, dv)
617 def _planform_variable(self, key):
618 if key in self.variables:
619 return self.variables[key]
620 return csdl.Variable(value=np.array([self.baseline[key]]))
622 def sectional(self, sectional):
623 sweep = self._planform_variable('sweep')
624 taper_ratio = self._planform_variable('taper_ratio')
625 if self.derived == 'span':
626 aspect_ratio = self._planform_variable('aspect_ratio')
627 root_chord = self._planform_variable('root_chord')
628 semi_span = aspect_ratio * root_chord * (1. + taper_ratio) / 4.
629 else:
630 semi_span = self._planform_variable('span') / 2.
631 if self.derived == 'root_chord':
632 aspect_ratio = self._planform_variable('aspect_ratio')
633 root_chord = 4. * semi_span / (aspect_ratio * (1. + taper_ratio))
634 else:
635 root_chord = self._planform_variable('root_chord')
636 aspect_ratio = 4. * semi_span / (root_chord * (1. + taper_ratio))
637 chord = root_chord * (1. - (1. - taper_ratio) * self.eta) # c_j'
638 ratio = chord / self.c0 # r_j
639 y_new = semi_span * self.eta
640 x_ref = self.baseline['x_ref_root'] + y_new * csdl.tan(sweep)
642 # chord: scale about the reference chord line, and about the chord
643 # plane (t/c kept)
644 sectional.add_stretch((ratio - 1.) * self.extent_x, 0, pivot=self.pivots)
645 sectional.add_stretch((ratio - 1.) * self.extent_z, 2, pivot=self.pivots)
646 sectional.add_translation(x_ref - self.x_ref0, [1., 0., 0.])
647 sectional.add_translation(y_new - self.y, [0., 1., 0.])
648 # thickness modes scale the (chord-scaled) height above the chord
649 # plane; camber modes move by a fraction of the new chord
650 n1 = self.parameterization.ffd.shape[1]
651 for key, scale in (('thickness', ratio), ('camber', chord)):
652 if key not in self.variables:
653 continue
654 modes = self.section_spec[key]['modes']
655 values = self.variables[key].reshape((n1, modes))
656 for m in range(modes):
657 sectional.add_mode(values[:, m] * scale, [0., 0., 1.],
658 self.mode_weights[key][m])
660 tip_chord = root_chord * taper_ratio
661 self._outputs = {'sweep': sweep, 'aspect_ratio': aspect_ratio,
662 'root_chord': root_chord, 'taper_ratio': taper_ratio,
663 'tip_chord': tip_chord, 'semi_span': semi_span,
664 'span': 2. * semi_span,
665 'area': (root_chord + tip_chord) * semi_span,
666 'section_chords': chord}
668 def outputs(self):
669 """Planform quantities as csdl variables: sweep, aspect_ratio,
670 root_chord, taper_ratio, tip_chord, semi_span, span, area (full wing)
671 and section_chords (at the FFD sections)."""
672 return getattr(self, '_outputs', {})
674 def describe(self):
675 b = self.baseline
676 lines = [" wing baseline: sweep({:.2f} c) {:.4f} deg, AR {:.5g}, root chord {:.5g}, "
677 "taper {:.5g}, span {:.5g}, area {:.5g} (edge fit residual {:.2e}); "
678 "{} derived".format(self.f, np.degrees(b['sweep']), b['aspect_ratio'],
679 b['root_chord'], b['taper_ratio'], b['span'], b['area'],
680 b['fit_residual'], self.derived)]
681 if self.section_spec:
682 n1 = self.parameterization.ffd.shape[1]
683 lines.append(" wing sections: {} spanwise stations at y = {}; {} (max chordwise "
684 "mode fit residual {:.2e})".format(
685 n1, np.array2string(self.y, precision=3),
686 ", ".join("{} {} modes".format(key, spec['modes'])
687 for key, spec in self.section_spec.items()),
688 self.mode_fit_residual))
689 return lines