Coverage for dragonfly_sim/core/shape_parameterization.py: 92%
353 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1"""Wall-shape parameterization built on lsdo_geo.
3Maps shape design variables to the motion of the mesh wall nodes, and measures
4the deformed wall for geometric constraints. Three pieces:
6* ``WallFFD`` -- an ``lsdo_geo.FFDBlock`` around the wall nodes. Its output,
7 ``wall_displacement(C)``, is the replicated ``(M_global, dim)`` array that
8 ``IDWarp_jax.evaluate`` consumes.
9* Shape layers that produce the block's coefficients ``C``:
10 ``SectionalShape`` (lsdo_geo ``SectionalParameterization``: camber/thickness
11 per chordwise station in 2D; twist/chord/sweep/dihedral/span per spanwise
12 section in 3D; plus additive non-affine sectional modes) and
13 ``WallFFD.apply_cp_motions`` (raw control-point motions). Compose them as
14 ``C = apply_cp_motions(SectionalShape.apply(C0), cp_motions)`` -- see
15 ``SectionalShape`` for why the sectional layer must act on the constant
16 baseline. ``shape_design`` packages these into design-variable layers.
17* ``WallGeometry`` -- geometric constraint quantities of the deformed wall:
18 enclosed area (2D) / volume (3D), thickness at stations, and 3D planform
19 (section chords, span, planform area).
21Everything here is csdl_alpha graph code on replicated arrays, so derivatives
22come from csdl's autodiff and every rank computes identical values.
23"""
24import contextlib
25import itertools
27import numpy as np
28import scipy.sparse as sps
29import csdl_alpha as csdl
30import lsdo_function_spaces as lfs
31import lsdo_geo as lg
33from dragonfly_sim.utils.mpi_utils import allgather_numpy
36@contextlib.contextmanager
37def _inline_recording():
38 """Record the enclosed csdl operations inline (values computed as built).
40 lsdo_geo reads intermediate ``.value``'s while building a sectional
41 parameterization (stretch axes and origins, integer-axis directions), which
42 a non-inline recorder -- the drivers' -- leaves as None. Every input of
43 such a subgraph has a value (constant baselines and design variables), so
44 the subgraph alone can be recorded inline.
45 """
46 recorder = csdl.get_current_recorder()
47 was_inline = recorder.inline
48 recorder.inline = True
49 try:
50 yield
51 finally:
52 recorder.inline = was_inline
55def validate_cp_coord_opt_idxs(cp_coord_opt_idxs, dim):
56 """Check and return the optimized control-point directions as an int array.
58 Strictly increasing is required, not merely unique: the cp_motions design
59 variable's trailing axis is a compressed index into this array (see
60 ``WallFFD.apply_cp_motions``), so a fixed order is what lets the bounds
61 built by ``ffd_dv_utils.build_cp_motion_dv`` line up with the DV columns
62 without either side having to know how the other was written.
63 """
64 if cp_coord_opt_idxs is None:
65 return np.arange(dim)
66 idxs = np.asarray(cp_coord_opt_idxs)
67 if not np.issubdtype(idxs.dtype, np.integer):
68 raise TypeError("cp_coord_opt_idxs must be integer spatial direction "
69 "indices, got dtype {}".format(idxs.dtype))
70 if idxs.ndim != 1:
71 raise ValueError("cp_coord_opt_idxs must be 1-D, got shape {}".format(idxs.shape))
72 if idxs.size == 0:
73 raise ValueError("cp_coord_opt_idxs is empty: at least one coordinate "
74 "direction must be a design variable")
75 if np.any(np.diff(idxs) <= 0):
76 raise ValueError("cp_coord_opt_idxs must be strictly increasing, got {}"
77 .format(idxs.tolist()))
78 if idxs[0] < 0 or idxs[-1] >= dim:
79 raise ValueError("cp_coord_opt_idxs must lie in [0, {}) for a {}-D mesh, got {}"
80 .format(dim, dim, idxs.tolist()))
81 return idxs
84def ffd_corners_from_bounding_box(coords, margin=1e-4):
85 """Axis-aligned FFD block corners around `coords`, widened by margin*extent.
87 Returns the ``(2,)*dim + (dim,)`` corner array that
88 ``lfs.create_b_spline_from_corners`` takes: ``corners[i0, i1, ..., d]`` is
89 the min (i_d = 0) or max (i_d = 1) of coordinate d.
90 """
91 coords = np.asarray(coords, dtype=np.float64)
92 dim = coords.shape[1]
93 lo = coords.min(axis=0)
94 hi = coords.max(axis=0)
95 delta = hi - lo
96 bounds = np.stack([lo - margin * delta, hi + margin * delta]) # (2, dim)
97 corners = np.zeros((2,) * dim + (dim,))
98 for idx in itertools.product((0, 1), repeat=dim):
99 corners[idx] = [bounds[idx[d], d] for d in range(dim)]
100 return corners
103def ffd_corners_from_corner_list(ffd_block_corner_list, margin=1e-4):
104 """FFD block corners lofted between two rectangular end surfaces.
106 ``ffd_block_corner_list`` holds exactly two surfaces, each given by two
107 opposite corners, e.g. the root and tip rectangles of a wing block. The
108 "connecting" direction is the axis along which the surface midpoints differ
109 most (spanwise for a wing); each surface spans the remaining axes. Each
110 surface is widened by margin*extent along the axes it spans. The resulting
111 block is tapered/swept, and ``lfs.create_b_spline_from_corners`` fills it
112 with the same trilinear lattice the pre-lsdo_geo construction produced.
113 """
114 if len(ffd_block_corner_list) != 2:
115 raise ValueError("Expected exactly 2 surface corner pairs, got {}"
116 .format(len(ffd_block_corner_list)))
117 dim = len(ffd_block_corner_list[0][0])
118 surface_mins = np.zeros((2, dim))
119 surface_maxs = np.zeros((2, dim))
120 for s, (corner_a, corner_b) in enumerate(ffd_block_corner_list):
121 corner_a, corner_b = np.asarray(corner_a, float), np.asarray(corner_b, float)
122 surface_mins[s] = np.minimum(corner_a, corner_b)
123 surface_maxs[s] = np.maximum(corner_a, corner_b)
125 midpoints = 0.5 * (surface_mins + surface_maxs)
126 connecting_dir = int(np.argmax(np.abs(midpoints[1] - midpoints[0])))
128 delta = surface_maxs - surface_mins
129 widen = np.where(delta > 0, margin * delta, 0.)
130 surface_mins = surface_mins - widen
131 surface_maxs = surface_maxs + widen
133 corners = np.zeros((2,) * dim + (dim,))
134 for idx in itertools.product((0, 1), repeat=dim):
135 s = idx[connecting_dir]
136 corners[idx] = [surface_mins[s, d] if (d == connecting_dir or idx[d] == 0)
137 else surface_maxs[s, d] for d in range(dim)]
138 return corners
141class WallFFD:
142 """lsdo_geo FFD block embedding the mesh wall nodes.
144 Replicated across ranks: every rank evaluates the FFD at all ``M_global``
145 wall nodes and returns the full global array, rather than evaluating only
146 the nodes it owns and having the pieces assembled afterwards.
148 That is deliberate. A rank owning no wall nodes would otherwise return a
149 ``(0, dim)`` array, and the reverse chain back to the design variables would
150 reduce an empty array -- which XLA constant-folds at compile time. The
151 downstream gradient reduction would then have no data dependence on
152 anything on that rank, and since ``jax.pure_callback`` is declared
153 effect-free, its collective could be scheduled arbitrarily early relative
154 to the other ranks'. That is an ordering deadlock the ranks cannot recover
155 from. Evaluating the full wall everywhere keeps every rank's chain
156 non-degenerate and lets the gradient reduction happen inside
157 ``IDWarp_jax.compute_jacvec_product``, on an array every rank genuinely
158 computes.
160 Only the projection is rank-local: each rank projects the wall nodes it
161 owns (``local_rows``) and the resulting parametric coordinates are gathered
162 in rank order.
164 Parameters
165 ----------
166 comm : mpi4py communicator
167 global_wall_coords : (M_global, dim) array
168 All wall node coordinates, in the rank-block order of
169 ``Mesh.boundary_node_set((wall_tag,)).coords``.
170 local_rows : int array
171 The contiguous rows of ``global_wall_coords`` this rank owns
172 (that node set's ``global_idx``).
173 shape : sequence of int
174 Control points per parametric direction.
175 degree : int or sequence of int
176 ffd_block_corner_list : optional
177 Two-surface block definition, see ``ffd_corners_from_corner_list``.
178 None wraps the wall's bounding box.
179 margin : float
180 Relative widening of the block beyond the wall / the given surfaces.
181 projection_newton_tol : float
182 Newton tolerance of the wall node projection into the block.
183 """
185 def __init__(self, comm, global_wall_coords, local_rows, shape, degree=2,
186 ffd_block_corner_list=None, margin=1e-4, projection_newton_tol=1e-12,
187 name='wall_ffd'):
188 self.comm = comm
189 self.baseline_wall_coords = np.ascontiguousarray(global_wall_coords, dtype=np.float64)
190 self.M_global, self.dim = self.baseline_wall_coords.shape
191 self.shape = tuple(int(n) for n in shape)
192 if len(self.shape) != self.dim:
193 raise ValueError("FFD shape {} does not match the {}-D wall"
194 .format(self.shape, self.dim))
195 degree = (degree,) * self.dim if np.isscalar(degree) else tuple(int(p) for p in degree)
196 self.degree = degree
198 if ffd_block_corner_list is None:
199 corners = ffd_corners_from_bounding_box(self.baseline_wall_coords, margin)
200 self.ffd_block_mode = 'automatic'
201 else:
202 corners = ffd_corners_from_corner_list(ffd_block_corner_list, margin)
203 self.ffd_block_mode = 'manual'
204 self.corners = corners
206 b_spline = lfs.create_b_spline_from_corners(
207 corners, degree=self.degree, num_coefficients=self.shape, name=name)
208 # No embedded_entities here: FFDBlock.embed_entities would project with
209 # refine_projection, which assumes 3-D points, and would pickle the
210 # result into ./stored_files/projections. The wall and the probe sets
211 # are projected below instead, and appended as embedded entities.
212 self.block = lg.FFDBlock(space=b_spline.space, coefficients=b_spline.coefficients,
213 name=name)
214 self.baseline_coefficients = np.array(self.block.coefficients.value, dtype=np.float64)
216 local_rows = np.asarray(local_rows, dtype=np.int64)
217 if local_rows.size and np.any(np.diff(local_rows) != 1):
218 raise ValueError("local_rows must be a contiguous block of wall rows")
219 self.projection_newton_tol = projection_newton_tol
220 local_u = self._project(self.baseline_wall_coords[local_rows])
221 self.wall_parametric_coordinates = allgather_numpy(
222 comm, local_u.ravel(), self.M_global * self.dim,
223 ).reshape(self.M_global, self.dim)
224 self.block.embedded_entity_parametric_coordinates.append(self.wall_parametric_coordinates)
226 self._probe_points = []
228 def _project(self, points):
229 if points.shape[0] == 0:
230 return np.zeros((0, self.dim))
231 return np.asarray(self.block.project(
232 points=points, plot=False, force_reproject=True,
233 grid_search_density_parameter=10, grid_search_subtraction_cutoff=1000,
234 newton_tolerance=self.projection_newton_tol, projection_tolerance=None,
235 do_pickles=False), dtype=np.float64).reshape(-1, self.dim)
237 def add_probe_points(self, points):
238 """Embed extra points (replicated input); returns their probe index.
240 Projected on rank 0 and broadcast, so every rank holds identical
241 parametric coordinates.
242 """
243 points = np.ascontiguousarray(points, dtype=np.float64).reshape(-1, self.dim)
244 u = self._project(points) if self.comm.rank == 0 else None
245 u = self.comm.bcast(u, root=0)
246 self.block.embedded_entity_parametric_coordinates.append(u)
247 self._probe_points.append((points, u))
248 return len(self._probe_points) - 1
250 def baseline_variable(self):
251 """The baseline coefficients as a constant csdl Variable."""
252 return csdl.Variable(value=self.baseline_coefficients.copy())
254 def apply_cp_motions(self, coefficients, cp_motions, cp_coord_opt_idxs):
255 """Add raw control-point motions to `coefficients`.
257 ``cp_motions`` has shape ``ffd_shape + (len(cp_coord_opt_idxs),)``;
258 column i moves spatial component ``cp_coord_opt_idxs[i]`` and the other
259 directions are left unchanged. The index tuples are built rather than
260 written as literals so this works for any dimension.
261 """
262 idxs = validate_cp_coord_opt_idxs(cp_coord_opt_idxs, self.dim)
263 expected_shape = self.shape + (idxs.size,)
264 if tuple(cp_motions.shape) != expected_shape:
265 raise ValueError(
266 "cp_motions has shape {}, expected {} (ffd_shape + (n_opt,), where "
267 "n_opt = len(cp_coord_opt_idxs) = {}). See ffd_dv_utils.build_cp_motion_dv "
268 "for constructing a design variable of the right shape."
269 .format(tuple(cp_motions.shape), expected_shape, idxs.size))
270 updated = coefficients
271 for i, coord in enumerate(idxs):
272 cp_index = tuple([slice(None)] * self.dim + [int(coord)])
273 motion_index = tuple([slice(None)] * self.dim + [i])
274 updated = updated.set(csdl.slice[cp_index],
275 coefficients[cp_index] + cp_motions[motion_index])
276 return updated
278 def _displacement_at(self, parametric_coordinates, coefficients):
279 # The FFD map is linear in the coefficients, so the displacement is the
280 # map applied to C - C0: exact zero at the baseline, whatever the
281 # projection residual of the baseline points.
282 delta = coefficients - self.baseline_variable()
283 return self.block.evaluate(parametric_coordinates, coefficients=delta)
285 def wall_displacement(self, coefficients):
286 """Replicated ``(M_global, dim)`` wall node displacement for `coefficients`."""
287 return self._displacement_at(self.wall_parametric_coordinates, coefficients)
289 def probe_positions(self, coefficients, probe):
290 """Deformed positions of probe set `probe` (index from add_probe_points)."""
291 points, u = self._probe_points[probe]
292 return self._displacement_at(u, coefficients) + points
295class SectionalShape:
296 """Sectional shape variables on a WallFFD, via lsdo_geo's SectionalParameterization.
298 The block is seen as a stack of sections normal to ``principal_dim``: in 2D
299 use principal_dim = 0 (chordwise stations, each a vertical column of control
300 points); for a wing, the spanwise parametric direction. Each variable holds
301 either one value per section, or -- when it has fewer entries -- the
302 coefficients of a B-spline profile over the sections.
304 The sectional layer must act on the *constant* baseline coefficients.
305 lsdo_geo takes stretch and translation axes, and stretch origins, from the
306 points' values when the graph is built, so feeding it DV-dependent points
307 would bake the initial DV values into the map. Add raw control-point
308 motions afterwards (``WallFFD.apply_cp_motions``).
310 Named presets (``add``), with lsdo_geo's conventions: a stretch is the
311 additive change in section extent along the axis (linear across the
312 section, zero at the pivot); a rotation is in radians about the axis
313 through the pivot; a translation is a length. A pivot is a parametric
314 coordinate within the section, either one for all sections or one row per
315 section.
317 lsdo_geo's operations are affine within each section. ``add_mode`` adds a
318 non-affine one -- a fixed control-point displacement pattern per section,
319 scaled by a sectional amplitude (e.g. a camber line) -- on top of the
320 lsdo_geo result; it is linear in the amplitude and reads no values.
322 ========== ======================================== =======
323 kind operation dims
324 ========== ======================================== =======
325 camber translation along the vertical axis 2D
326 thickness stretch along the vertical param. axis 2D, 3D
327 twist rotation about the principal axis 3D
328 chord stretch along parametric axis 0 3D
329 sweep translation along x 3D
330 dihedral translation along z 3D
331 span translation along y 3D
332 ========== ======================================== =======
333 """
335 def __init__(self, ffd, principal_dim=0):
336 self.ffd = ffd
337 self.dim = ffd.dim
338 if not 0 <= principal_dim < self.dim:
339 raise ValueError("principal_dim must lie in [0, {})".format(self.dim))
340 self.principal_dim = principal_dim
341 self.num_sections = ffd.shape[principal_dim]
342 # (kind, values, axis, pivot, profile_degree); turned into lsdo_geo
343 # SectionalParameters inside apply(), see there for why
344 self._specs = []
345 # (values, (n_coefficients, num_sections) map, profile_degree)
346 self._modes = []
348 def _per_section(self, values, profile_degree=2):
349 n = int(np.prod(values.shape))
350 if n == self.num_sections:
351 return values.reshape((self.num_sections,))
352 if not 1 <= n < self.num_sections:
353 raise ValueError("a sectional variable needs {} entries (one per section) "
354 "or fewer (profile coefficients), got {}"
355 .format(self.num_sections, n))
356 space = lfs.BSplineSpace(num_parametric_dimensions=1,
357 degree=min(profile_degree, n - 1), coefficients_shape=(n,))
358 basis = space.compute_basis_matrix(
359 np.linspace(0., 1., self.num_sections).reshape(-1, 1))
360 basis = basis.toarray() if sps.issparse(basis) else np.asarray(basis)
361 return csdl.matvec(basis, values.reshape((n,)))
363 def _section_pivot(self, pivot):
364 """One pivot for all sections, (n_section_dims,), or one per section,
365 (num_sections, n_section_dims)."""
366 n_section_dims = self.dim - 1
367 if pivot is None:
368 return None
369 pivot = np.asarray(pivot, dtype=np.float64)
370 if pivot.ndim == 2:
371 if pivot.shape != (self.num_sections, n_section_dims):
372 raise ValueError("per-section pivots must have shape {}, got {}".format(
373 (self.num_sections, n_section_dims), pivot.shape))
374 return pivot
375 pivot = pivot.reshape(-1)
376 if pivot.size != n_section_dims:
377 raise ValueError("pivot is a parametric coordinate within the section "
378 "({} entries), got {}".format(n_section_dims, pivot.size))
379 return pivot
381 def _check_size(self, values):
382 n = int(np.prod(values.shape))
383 if not 1 <= n <= self.num_sections:
384 raise ValueError("a sectional variable needs {} entries (one per section) "
385 "or fewer (profile coefficients), got {}"
386 .format(self.num_sections, n))
388 def add_translation(self, values, direction, profile_degree=2):
389 self._check_size(values)
390 self._specs.append(('translation', values, np.asarray(direction, dtype=np.float64),
391 None, profile_degree))
393 def add_stretch(self, values, parametric_axis, pivot=None, profile_degree=2):
394 self._check_size(values)
395 self._specs.append(('stretch', values, int(parametric_axis),
396 self._section_pivot(pivot), profile_degree))
398 def add_rotation(self, values, parametric_axis=None, pivot=None, profile_degree=2):
399 if self.dim != 3:
400 raise ValueError("sectional rotations need a 3-D block (lsdo_geo "
401 "rotates with quaternions)")
402 self._check_size(values)
403 axis = self.principal_dim if parametric_axis is None else int(parametric_axis)
404 self._specs.append(('rotation', values, axis, self._section_pivot(pivot), profile_degree))
406 def add_mode(self, values, direction, weights, profile_degree=2):
407 """Add ``values[j] * weights[..., j-th section] * direction`` to the
408 control points.
410 weights has the block's shape (``ffd.shape``); the entries of section j
411 (index j along principal_dim) are that section's displacement pattern
412 per unit amplitude. direction is a physical vector.
413 """
414 self._check_size(values)
415 weights = np.asarray(weights, dtype=np.float64)
416 if weights.shape != self.ffd.shape:
417 raise ValueError("mode weights must have the block's shape {}, got {}"
418 .format(self.ffd.shape, weights.shape))
419 direction = np.asarray(direction, dtype=np.float64).reshape(-1)
420 if direction.size != self.dim:
421 raise ValueError("mode direction must have {} entries".format(self.dim))
422 section_index = np.indices(self.ffd.shape)[self.principal_dim]
423 mode_map = np.zeros(self.ffd.shape + (self.dim, self.num_sections))
424 for j in range(self.num_sections):
425 mode_map[..., j] = np.where(section_index == j, weights, 0.)[..., None] * direction
426 self._modes.append((values, mode_map.reshape(-1, self.num_sections), profile_degree))
428 def add(self, kind, values, pivot=None, profile_degree=2):
429 """Add a named sectional variable; see the class docstring."""
430 vertical = self.dim - 1
431 unit = np.eye(self.dim)
432 if kind == 'camber' and self.dim == 2:
433 self.add_translation(values, unit[vertical], profile_degree)
434 elif kind == 'thickness':
435 self.add_stretch(values, vertical, pivot, profile_degree)
436 elif kind == 'twist':
437 self.add_rotation(values, None, pivot, profile_degree)
438 elif kind == 'chord' and self.dim == 3:
439 self.add_stretch(values, 0, pivot, profile_degree)
440 elif kind == 'sweep' and self.dim == 3:
441 self.add_translation(values, unit[0], profile_degree)
442 elif kind == 'dihedral' and self.dim == 3:
443 self.add_translation(values, unit[2], profile_degree)
444 elif kind == 'span' and self.dim == 3:
445 self.add_translation(values, unit[1], profile_degree)
446 else:
447 raise ValueError("unknown sectional variable {!r} for a {}-D block"
448 .format(kind, self.dim))
450 def apply(self, coefficients):
451 """Sectionally deformed coefficients. Pass ``ffd.baseline_variable()``."""
452 if not self._specs and not self._modes:
453 return coefficients
454 # lsdo_geo reads values while it builds the subgraph; record it, profiles
455 # included, inline (see _inline_recording).
456 with _inline_recording():
457 updated = coefficients
458 if self._specs:
459 parameters = lg.SectionalParameters()
460 for kind, values, axis, pivot, profile_degree in self._specs:
461 per_section = self._per_section(values, profile_degree)
462 if kind == 'translation':
463 parameters.add_translation(axis, per_section)
464 elif kind == 'stretch':
465 parameters.add_stretch(axis, per_section, pivot)
466 else:
467 parameters.add_rotation(axis, per_section, pivot)
468 parameterization = lg.SectionalParameterization(
469 parameterized_points=coefficients,
470 principal_parametric_dimension=self.principal_dim)
471 updated = parameterization.evaluate(parameters)
472 for values, mode_map, profile_degree in self._modes:
473 per_section = self._per_section(values, profile_degree)
474 updated = updated + csdl.matvec(mode_map, per_section).reshape(updated.shape)
475 return updated
478def _selection_matrix(ids, n_cols):
479 ids = np.asarray(ids, dtype=np.int64)
480 return sps.csr_matrix((np.ones(ids.size), (np.arange(ids.size), ids)),
481 shape=(ids.size, n_cols))
484class WallGeometry:
485 """Geometric constraint quantities of the FFD-deformed wall.
487 Parameters
488 ----------
489 ffd : WallFFD
490 wall_facets : (n_facets, k) int array
491 Replicated wall facets as rows of wall node ids in ``[0, M_global)``
492 (``Mesh.wall_facets(wall_tag)``): segments in 2D, triangles or quads
493 (cyclic vertex order) in 3D. Their winding only needs to be consistent;
494 the sign of the enclosed measure is fixed from the baseline.
496 In 3D the wall may be open on a symmetry plane through the origin (y = 0):
497 that plane contributes nothing to the divergence-theorem volume, as long as
498 the root section keeps y = 0.
499 """
501 def __init__(self, ffd, wall_facets):
502 self.ffd = ffd
503 self.dim = ffd.dim
504 facets = np.asarray(wall_facets, dtype=np.int64)
505 if facets.ndim != 2 or facets.size == 0:
506 raise ValueError("wall_facets must be a non-empty (n_facets, k) array")
507 if facets.min() < 0 or facets.max() >= ffd.M_global:
508 raise ValueError("wall_facets must index wall nodes in [0, {})".format(ffd.M_global))
509 if self.dim == 2:
510 if facets.shape[1] != 2:
511 raise ValueError("2D wall facets must be segments")
512 self.simplices = facets
513 else:
514 if facets.shape[1] == 3:
515 self.simplices = facets
516 elif facets.shape[1] == 4:
517 self.simplices = np.concatenate([facets[:, [0, 1, 2]], facets[:, [0, 2, 3]]])
518 else:
519 raise ValueError("3D wall facets must be triangles or quads")
520 self._selectors = [_selection_matrix(self.simplices[:, j], ffd.M_global)
521 for j in range(self.simplices.shape[1])]
522 raw = self._raw_measure_numpy(ffd.baseline_wall_coords)
523 if not abs(raw) > 0.:
524 raise ValueError("the baseline wall encloses no area/volume; "
525 "are the wall facets complete?")
526 self._sign = float(np.sign(raw))
527 self.baseline_measure = abs(raw)
529 # ---- enclosed area / volume -------------------------------------------
530 def _raw_measure_numpy(self, X):
531 v = [X[self.simplices[:, j]] for j in range(self.simplices.shape[1])]
532 if self.dim == 2:
533 return 0.5 * np.sum(v[0][:, 0] * v[1][:, 1] - v[1][:, 0] * v[0][:, 1])
534 return np.sum(np.einsum('ij,ij->i', v[0], np.cross(v[1], v[2]))) / 6.
536 def _gather(self, X, j):
537 cols = [csdl.sparse.matvec(self._selectors[j], X[:, d].reshape((self.ffd.M_global, 1)))
538 for d in range(self.dim)]
539 return csdl.concatenate(cols, axis=1)
541 def enclosed_measure(self, wall_displacement):
542 """Enclosed area (2D) or volume (3D) of the deformed wall, positive."""
543 X = wall_displacement + self.ffd.baseline_wall_coords
544 v = [self._gather(X, j) for j in range(self.simplices.shape[1])]
545 if self.dim == 2:
546 raw = 0.5 * csdl.sum(v[0][:, 0] * v[1][:, 1] - v[1][:, 0] * v[0][:, 1])
547 else:
548 raw = csdl.sum(v[0] * csdl.cross(v[1], v[2], axis=1)) / 6.
549 return self._sign * raw
551 # ---- baseline-geometry queries (numpy, replicated) ---------------------
552 def _vertical_line_hits(self, x_in_plane):
553 """Vertical-coordinate values where the vertical line through
554 `x_in_plane` (x in 2D, (x, y) in 3D) crosses the baseline wall."""
555 X = self.ffd.baseline_wall_coords
556 if self.dim == 2:
557 a, b = X[self.simplices[:, 0]], X[self.simplices[:, 1]]
558 x0 = float(x_in_plane)
559 crosses = (a[:, 0] - x0) * (b[:, 0] - x0) <= 0.
560 crosses &= a[:, 0] != b[:, 0]
561 t = (x0 - a[crosses, 0]) / (b[crosses, 0] - a[crosses, 0])
562 return a[crosses, 1] + t * (b[crosses, 1] - a[crosses, 1])
563 p = np.asarray(x_in_plane, dtype=np.float64)
564 a, b, c = (X[self.simplices[:, j]] for j in range(3))
565 # barycentric coordinates in the xy projection
566 v0, v1, v2 = b[:, :2] - a[:, :2], c[:, :2] - a[:, :2], p - a[:, :2]
567 den = v0[:, 0] * v1[:, 1] - v1[:, 0] * v0[:, 1]
568 ok = np.abs(den) > 1e-14 * (np.abs(v0).max() * np.abs(v1).max() + 1e-300)
569 den = np.where(ok, den, 1.)
570 l1 = (v2[:, 0] * v1[:, 1] - v1[:, 0] * v2[:, 1]) / den
571 l2 = (v0[:, 0] * v2[:, 1] - v2[:, 0] * v0[:, 1]) / den
572 tol = 1e-12
573 inside = ok & (l1 >= -tol) & (l2 >= -tol) & (l1 + l2 <= 1. + tol)
574 return (a[inside, 2] + l1[inside] * (b[inside, 2] - a[inside, 2])
575 + l2[inside] * (c[inside, 2] - a[inside, 2]))
577 def _section_leading_trailing_edges(self, y0):
578 """Baseline (LE, TE) points of the 3D wall's section at span y0: the
579 min-x and max-x points where the plane y = y0 cuts the wall."""
580 X = self.ffd.baseline_wall_coords
581 pts = []
582 for j, k in ((0, 1), (1, 2), (2, 0)):
583 a, b = X[self.simplices[:, j]], X[self.simplices[:, k]]
584 crosses = ((a[:, 1] - y0) * (b[:, 1] - y0) <= 0.) & (a[:, 1] != b[:, 1])
585 t = (y0 - a[crosses, 1]) / (b[crosses, 1] - a[crosses, 1])
586 pts.append(a[crosses] + t[:, None] * (b[crosses] - a[crosses]))
587 pts = np.concatenate(pts)
588 if pts.shape[0] == 0:
589 raise ValueError("the plane y = {} does not cut the wall".format(y0))
590 return pts[np.argmin(pts[:, 0])], pts[np.argmax(pts[:, 0])]
592 def section_edges(self, y0):
593 """Baseline leading- and trailing-edge points of the 3D wall's section
594 at span y0 (the min-x and max-x points of the cut y = y0)."""
595 if self.dim != 3:
596 raise ValueError("section_edges is for 3D walls")
597 return self._section_leading_trailing_edges(y0)
599 def _chord_range(self, y0=None):
600 if self.dim == 2:
601 x = self.ffd.baseline_wall_coords[:, 0]
602 return x.min(), x.max()
603 le, te = self._section_leading_trailing_edges(y0)
604 return le[0], te[0]
606 # ---- thickness at stations --------------------------------------------
607 def add_thickness_stations(self, chord_fractions, span_stations=None):
608 """Embed upper/lower wall points at chordwise stations; returns a handle.
610 2D: one station per chord fraction. 3D: the tensor product of chord
611 fractions (of the local section chord) and span stations (y values).
612 """
613 fractions = np.atleast_1d(np.asarray(chord_fractions, dtype=np.float64))
614 upper, lower = [], []
615 spans = [None] if self.dim == 2 else np.atleast_1d(span_stations)
616 for y0 in spans:
617 x_le, x_te = self._chord_range(y0)
618 for f in fractions:
619 x0 = x_le + f * (x_te - x_le)
620 where = x0 if self.dim == 2 else np.array([x0, y0])
621 hits = self._vertical_line_hits(where)
622 if hits.size < 2:
623 raise ValueError("thickness station at {} does not cross the wall "
624 "twice".format(where))
625 base = [x0] if self.dim == 2 else [x0, y0]
626 upper.append(base + [hits.max()])
627 lower.append(base + [hits.min()])
628 upper, lower = np.array(upper), np.array(lower)
629 handle = {'upper': self.ffd.add_probe_points(upper),
630 'lower': self.ffd.add_probe_points(lower),
631 'baseline': np.linalg.norm(upper - lower, axis=1)}
632 return handle
634 def thickness(self, coefficients, handle):
635 """Deformed thickness at the stations of `handle`."""
636 up = self.ffd.probe_positions(coefficients, handle['upper'])
637 lo = self.ffd.probe_positions(coefficients, handle['lower'])
638 return csdl.norm(up - lo, axes=(1,))
640 def thickness_ratio(self, coefficients, handle):
641 """Deformed over baseline thickness at the stations of `handle`."""
642 return self.thickness(coefficients, handle) / handle['baseline']
644 # ---- 3D planform ------------------------------------------------------
645 def add_planform_stations(self, span_stations):
646 """Embed the LE and TE points of the wall sections at the given span
647 stations (increasing y, root first); returns a handle.
649 ``handle['baseline']`` holds the undeformed chord, span and area as
650 numbers, for constraint bounds (the csdl values are not available while
651 the graph is recorded).
652 """
653 if self.dim != 3:
654 raise ValueError("planform quantities are for 3D walls")
655 span_stations = np.asarray(span_stations, dtype=np.float64)
656 if span_stations.size < 2 or np.any(np.diff(span_stations) <= 0):
657 raise ValueError("need at least two increasing span stations")
658 edges = [self._section_leading_trailing_edges(y0) for y0 in span_stations]
659 le = np.array([e[0] for e in edges])
660 te = np.array([e[1] for e in edges])
661 chord = np.linalg.norm((te - le)[:, :2], axis=1)
662 y = le[:, 1]
663 baseline = {'chord': chord, 'span': float(y[-1] - y[0]),
664 'area': float(np.sum(0.5 * (chord[1:] + chord[:-1]) * np.diff(y)))}
665 return {'le': self.ffd.add_probe_points(le), 'te': self.ffd.add_probe_points(te),
666 'n': span_stations.size, 'baseline': baseline}
668 def planform(self, coefficients, handle):
669 """Deformed section chords, span and trapezoidal planform area.
671 Chords are LE-TE distances projected onto the xy plane; span and area
672 use the y coordinates of the leading edges.
673 """
674 le = self.ffd.probe_positions(coefficients, handle['le'])
675 te = self.ffd.probe_positions(coefficients, handle['te'])
676 chord = csdl.norm((te - le)[:, :2], axes=(1,))
677 y = le[:, 1]
678 n = handle['n']
679 dy = y[1:] - y[:n - 1]
680 area = csdl.sum(0.5 * (chord[1:] + chord[:n - 1]) * dy)
681 span = y[n - 1] - y[0]
682 return {'chord': chord, 'span': span, 'area': area}