Additional constraints (#975)

* Started assy solver refactoring

* First working version

* Added validation

+initial work on FixedPoint

* Fixed tests

* Unary constraints support and finish Fixed Point

* Added Fixed and FixedAxis

* Changed the locking logic

* Test for unary constraints

* Update cadquery/occ_impl/solver.py

Co-authored-by: Marcus Boyd <mwb@geosol.com.au>

* Simple validation test

* PointOnLine constraint

* PointOnLine test

* FixedRotation and FixedRotationAxis

* FixedRotation test

* Describe PointOnLine

* Add new constraints to the docs

* Use gp_Intrinsic_XYZ

* Apply suggestions from code review

Co-authored-by: Marcus Boyd <mwb@geosol.com.au>

* Apply suggestions from code review

Co-authored-by: Lorenz <hello@lorenz.space>

* Black fix

* Apply suggestions from code review

Co-authored-by: Lorenz <hello@lorenz.space>
Co-authored-by: Jeremy Wright <wrightjmf@gmail.com>

* Use windows-latest

Co-authored-by: Marcus Boyd <mwb@geosol.com.au>
Co-authored-by: Lorenz <hello@lorenz.space>
Co-authored-by: Jeremy Wright <wrightjmf@gmail.com>
This commit is contained in:
AU
2022-02-15 08:42:54 +01:00
committed by GitHub
parent ed8e62c783
commit 4108c7dc5f
5 changed files with 722 additions and 260 deletions

View File

@ -43,7 +43,7 @@ jobs:
- template: conda-build.yml@templates
parameters:
name: Windows
vmImage: 'vs2017-win2016'
vmImage: 'windows-latest'
py_maj: 3
py_min: ${{minor}}
conda_bld: 3.21.6

View File

@ -1,16 +1,18 @@
from functools import reduce
from typing import Union, Optional, List, Dict, Any, overload, Tuple, Iterator, cast
from typing_extensions import Literal
from typish import instance_of
from uuid import uuid1 as uuid
from .cq import Workplane
from .occ_impl.shapes import Shape, Compound, Face, Edge, Wire
from .occ_impl.geom import Location, Vector, Plane
from .occ_impl.shapes import Shape, Compound
from .occ_impl.geom import Location
from .occ_impl.assembly import Color
from .occ_impl.solver import (
ConstraintSolver,
ConstraintMarker,
Constraint as ConstraintPOD,
ConstraintSpec as Constraint,
UnaryConstraintKind,
BinaryConstraintKind,
)
from .occ_impl.exporters.assembly import (
exportAssembly,
@ -21,9 +23,6 @@ from .occ_impl.exporters.assembly import (
)
from .selectors import _expression_grammar as _selector_grammar
from OCP.BRepTools import BRepTools
from OCP.gp import gp_Pln, gp_Pnt
from OCP.Precision import Precision
# type definitions
AssemblyObjects = Union[Shape, Workplane, None]
@ -67,112 +66,6 @@ def _define_grammar():
_grammar = _define_grammar()
class Constraint(object):
"""
Geometrical constraint between two shapes of an assembly.
"""
objects: Tuple[str, ...]
args: Tuple[Shape, ...]
sublocs: Tuple[Location, ...]
kind: ConstraintKinds
param: Any
def __init__(
self,
objects: Tuple[str, ...],
args: Tuple[Shape, ...],
sublocs: Tuple[Location, ...],
kind: ConstraintKinds,
param: Any = None,
):
"""
Construct a constraint.
:param objects: object names referenced in the constraint
:param args: subshapes (e.g. faces or edges) of the objects
:param sublocs: locations of the objects (only relevant if the objects are nested in a sub-assembly)
:param kind: constraint kind
:param param: optional arbitrary parameter passed to the solver
"""
self.objects = objects
self.args = args
self.sublocs = sublocs
self.kind = kind
self.param = param
def _getAxis(self, arg: Shape) -> Vector:
if isinstance(arg, Face):
rv = arg.normalAt()
elif isinstance(arg, Edge) and arg.geomType() != "CIRCLE":
rv = arg.tangentAt()
elif isinstance(arg, Edge) and arg.geomType() == "CIRCLE":
rv = arg.normal()
else:
raise ValueError(f"Cannot construct Axis for {arg}")
return rv
def _getPln(self, arg: Shape) -> gp_Pln:
if isinstance(arg, Face):
rv = gp_Pln(self._getPnt(arg), arg.normalAt().toDir())
elif isinstance(arg, (Edge, Wire)):
normal = arg.normal()
origin = arg.Center()
plane = Plane(origin, normal=normal)
rv = plane.toPln()
else:
raise ValueError(f"Can not construct a plane for {arg}.")
return rv
def _getPnt(self, arg: Shape) -> gp_Pnt:
# check for infinite face
if isinstance(arg, Face) and any(
Precision.IsInfinite_s(x) for x in BRepTools.UVBounds_s(arg.wrapped)
):
# fall back to gp_Pln center
pln = arg.toPln()
center = Vector(pln.Location())
else:
center = arg.Center()
return center.toPnt()
def toPOD(self) -> ConstraintPOD:
"""
Convert the constraint to a representation used by the solver.
"""
rv: List[Tuple[ConstraintMarker, ...]] = []
for idx, (arg, loc) in enumerate(zip(self.args, self.sublocs)):
arg = arg.located(loc * arg.location())
if self.kind == "Axis":
rv.append((self._getAxis(arg).toDir(),))
elif self.kind == "Point":
rv.append((self._getPnt(arg),))
elif self.kind == "Plane":
rv.append((self._getAxis(arg).toDir(), self._getPnt(arg)))
elif self.kind == "PointInPlane":
if idx == 0:
rv.append((self._getPnt(arg),))
else:
rv.append((self._getPln(arg),))
else:
raise ValueError(f"Unknown constraint kind {self.kind}")
rv.append(self.param)
return cast(ConstraintPOD, tuple(rv))
class Assembly(object):
"""Nested assembly of Workplane and Shape objects defining their relative positions."""
@ -389,6 +282,12 @@ class Assembly(object):
) -> "Assembly":
...
@overload
def constrain(
self, q1: str, kind: ConstraintKinds, param: Any = None
) -> "Assembly":
...
@overload
def constrain(
self,
@ -401,12 +300,25 @@ class Assembly(object):
) -> "Assembly":
...
@overload
def constrain(
self, id1: str, s1: Shape, kind: ConstraintKinds, param: Any = None,
) -> "Assembly":
...
def constrain(self, *args, param=None):
"""
Define a new constraint.
"""
if len(args) == 3:
# dispatch on arguments
if len(args) == 2:
q1, kind = args
id1, s1 = self._query(q1)
elif len(args) == 3 and instance_of(args[1], UnaryConstraintKind):
q1, kind, param = args
id1, s1 = self._query(q1)
elif len(args) == 3:
q1, q2, kind = args
id1, s1 = self._query(q1)
id2, s2 = self._query(q2)
@ -421,11 +333,18 @@ class Assembly(object):
else:
raise ValueError(f"Incompatible arguments: {args}")
loc1, id1_top = self._subloc(id1)
loc2, id2_top = self._subloc(id2)
self.constraints.append(
Constraint((id1_top, id2_top), (s1, s2), (loc1, loc2), kind, param)
)
# handle unary and binary constraints
if instance_of(kind, UnaryConstraintKind):
loc1, id1_top = self._subloc(id1)
c = Constraint((id1_top,), (s1,), (loc1,), kind, param)
elif instance_of(kind, BinaryConstraintKind):
loc1, id1_top = self._subloc(id1)
loc2, id2_top = self._subloc(id2)
c = Constraint((id1_top, id2_top), (s1, s2), (loc1, loc2), kind, param)
else:
raise ValueError(f"Unknown constraint: {kind}")
self.constraints.append(c)
return self
@ -434,32 +353,53 @@ class Assembly(object):
Solve the constraints.
"""
# get all entities and number them
# Get all entities and number them. First entity is marked as locked
ents = {}
i = 0
lock_ix = 0
locked = []
for c in self.constraints:
for name in c.objects:
if name not in ents:
ents[name] = i
if name == self.name:
lock_ix = i
i += 1
if c.kind == "Fixed" or name == self.name:
locked.append(ents[name])
# Lock the first occuring entity if needed.
if not locked:
unary_objects = [
c.objects[0]
for c in self.constraints
if instance_of(c.kind, UnaryConstraintKind)
]
binary_objects = [
c.objects[0]
for c in self.constraints
if instance_of(c.kind, BinaryConstraintKind)
]
for b in binary_objects:
if b not in unary_objects:
locked.append(ents[b])
break
locs = [self.objects[n].loc for n in ents]
# construct the constraint mapping
constraints = []
for c in self.constraints:
constraints.append(((ents[c.objects[0]], ents[c.objects[1]]), c.toPOD()))
ixs = tuple(ents[obj] for obj in c.objects)
pods = c.toPODs()
for pod in pods:
constraints.append((ixs, pod))
# check if any constraints were specified
if not constraints:
raise ValueError("At least one constraint required")
# instantiate the solver
solver = ConstraintSolver(locs, constraints, locked=[lock_ix])
solver = ConstraintSolver(locs, constraints, locked=locked)
# solve
locs_new, self._solve_result = solver.solve()

View File

@ -1,17 +1,99 @@
from typing import Tuple, Union, Any, Callable, List, Optional, Dict
from typing import (
List,
Tuple,
Union,
Any,
Callable,
Optional,
Dict,
Literal,
cast as tcast,
Type,
)
from nptyping import NDArray as Array
from numpy import array, eye, zeros, pi
from math import radians
from typish import instance_of, get_type
from numpy import array, eye, pi
import nlopt
from OCP.gp import gp_Vec, gp_Pln, gp_Dir, gp_Pnt, gp_Trsf, gp_Quaternion
from OCP.gp import (
gp_Vec,
gp_Pln,
gp_Dir,
gp_Pnt,
gp_Trsf,
gp_Quaternion,
gp_XYZ,
gp_Lin,
gp_Intrinsic_XYZ,
)
from OCP.BRepTools import BRepTools
from OCP.Precision import Precision
from .geom import Location
from .geom import Location, Vector, Plane
from .shapes import Shape, Face, Edge, Wire
from ..types import Real
# type definitions
NoneType = type(None)
DOF6 = Tuple[float, float, float, float, float, float]
ConstraintMarker = Union[gp_Pln, gp_Dir, gp_Pnt]
ConstraintMarker = Union[gp_Pln, gp_Dir, gp_Pnt, gp_Lin, None]
UnaryConstraintKind = Literal[
"Fixed", "FixedPoint", "FixedAxis", "FixedRotation", "FixedRotationAxis"
]
BinaryConstraintKind = Literal["Plane", "Point", "Axis", "PointInPlane", "PointOnLine"]
ConstraintKind = Literal[
"Plane",
"Point",
"Axis",
"PointInPlane",
"Fixed",
"FixedPoint",
"FixedAxis",
"PointOnLine",
"FixedRotation",
"FixedRotationAxis",
]
# (arity, marker types, param type, conversion func)
ConstraintInvariants = {
"Point": (2, (gp_Pnt, gp_Pnt), Real, None),
"Axis": (
2,
(gp_Dir, gp_Dir),
Real,
lambda x: radians(x) if x is not None else None,
),
"PointInPlane": (2, (gp_Pnt, gp_Pln), Real, None),
"PointOnLine": (2, (gp_Pnt, gp_Lin), Real, None),
"Fixed": (1, (None,), Type[None], None),
"FixedPoint": (1, (gp_Pnt,), Tuple[Real, Real, Real], None),
"FixedAxis": (1, (gp_Dir,), Tuple[Real, Real, Real], None),
"FixedRotationAxis": (
1,
(None,),
Tuple[int, Real],
lambda x: (x[0], radians(x[1])),
),
}
# translation table for compound constraints {name : (name, ...), converter}
CompoundConstraints: Dict[
ConstraintKind, Tuple[Tuple[ConstraintKind, ...], Callable[[Any], Tuple[Any, ...]]]
] = {
"Plane": (("Axis", "Point"), lambda x: (radians(x) if x is not None else None, 0)),
"FixedRotation": (
("FixedRotationAxis", "FixedRotationAxis", "FixedRotationAxis"),
lambda x: tuple(enumerate(map(radians, x))),
),
}
# constraint POD type
Constraint = Tuple[
Tuple[ConstraintMarker, ...], Tuple[Optional[ConstraintMarker], ...], Optional[Any]
Tuple[ConstraintMarker, ...], ConstraintKind, Optional[Any],
]
NDOF = 6
@ -20,11 +102,299 @@ DIFF_EPS = 1e-10
TOL = 1e-12
MAXITER = 2000
# high-level constraint class - to be used by clients
class ConstraintSpec(object):
"""
Geometrical constraint specification between two shapes of an assembly.
"""
objects: Tuple[str, ...]
args: Tuple[Shape, ...]
sublocs: Tuple[Location, ...]
kind: ConstraintKind
param: Any
def __init__(
self,
objects: Tuple[str, ...],
args: Tuple[Shape, ...],
sublocs: Tuple[Location, ...],
kind: ConstraintKind,
param: Any = None,
):
"""
Construct a constraint.
:param objects: object names referenced in the constraint
:param args: subshapes (e.g. faces or edges) of the objects
:param sublocs: locations of the objects (only relevant if the objects are nested in a sub-assembly)
:param kind: constraint kind
:param param: optional arbitrary parameter passed to the solver
"""
# validate
if not instance_of(kind, ConstraintKind):
raise ValueError(f"Unknown constraint {kind}.")
if kind in CompoundConstraints:
kinds, convert_compound = CompoundConstraints[kind]
for k, p in zip(kinds, convert_compound(param)):
self._validate(args, k, p)
else:
self._validate(args, kind, param)
# convert here for simple constraints
convert = ConstraintInvariants[kind][-1]
param = convert(param) if convert else param
# store
self.objects = objects
self.args = args
self.sublocs = sublocs
self.kind = kind
self.param = param
def _validate(self, args: Tuple[Shape, ...], kind: ConstraintKind, param: Any):
arity, marker_types, param_type, converter = ConstraintInvariants[kind]
# check arity
if arity != len(args):
raise ValueError(
f"Invalid number of entities for constraint {kind}. Provided {len(args)}, required {arity}."
)
# check arguments
arg_check: Dict[Any, Callable[[Shape], Any]] = {
gp_Pnt: self._getPnt,
gp_Dir: self._getAxis,
gp_Pln: self._getPln,
gp_Lin: self._getLin,
None: lambda x: True, # dummy check for None marker
}
for a, t in zip(args, tcast(Tuple[Type[ConstraintMarker], ...], marker_types)):
try:
arg_check[t](a)
except ValueError:
raise ValueError(f"Unsupported entity {a} for constraint {kind}.")
# check parameter
if not instance_of(param, param_type) and param is not None:
raise ValueError(
f"Unsupported argument types {get_type(param)}, required {param_type}."
)
# check parameter conversion
try:
if param is not None and converter:
converter(param)
except Exception as e:
raise ValueError(f"Exception {e} occured in the parameter conversion")
def _getAxis(self, arg: Shape) -> gp_Dir:
if isinstance(arg, Face):
rv = arg.normalAt()
elif isinstance(arg, Edge) and arg.geomType() != "CIRCLE":
rv = arg.tangentAt()
elif isinstance(arg, Edge) and arg.geomType() == "CIRCLE":
rv = arg.normal()
else:
raise ValueError(f"Cannot construct Axis for {arg}")
return rv.toDir()
def _getPln(self, arg: Shape) -> gp_Pln:
if isinstance(arg, Face):
rv = gp_Pln(self._getPnt(arg), arg.normalAt().toDir())
elif isinstance(arg, (Edge, Wire)):
normal = arg.normal()
origin = arg.Center()
plane = Plane(origin, normal=normal)
rv = plane.toPln()
else:
raise ValueError(f"Cannot construct a plane for {arg}.")
return rv
def _getPnt(self, arg: Shape) -> gp_Pnt:
# check for infinite face
if isinstance(arg, Face) and any(
Precision.IsInfinite_s(x) for x in BRepTools.UVBounds_s(arg.wrapped)
):
# fall back to gp_Pln center
pln = arg.toPln()
center = Vector(pln.Location())
else:
center = arg.Center()
return center.toPnt()
def _getLin(self, arg: Shape) -> gp_Lin:
if isinstance(arg, (Edge, Wire)):
center = arg.Center()
tangent = arg.tangentAt()
else:
raise ValueError(f"Cannot construct a plane for {arg}.")
return gp_Lin(center.toPnt(), tangent.toDir())
def toPODs(self) -> Tuple[Constraint, ...]:
"""
Convert the constraint to a representation used by the solver.
NB: Compound constraints are decomposed into simple ones.
"""
# apply sublocation
args = tuple(
arg.located(loc * arg.location())
for arg, loc in zip(self.args, self.sublocs)
)
markers: List[Tuple[ConstraintMarker, ...]]
# convert to marker objects
if self.kind == "Axis":
markers = [(self._getAxis(args[0]), self._getAxis(args[1]),)]
elif self.kind == "Point":
markers = [(self._getPnt(args[0]), self._getPnt(args[1]))]
elif self.kind == "Plane":
markers = [
(self._getAxis(args[0]), self._getAxis(args[1]),),
(self._getPnt(args[0]), self._getPnt(args[1])),
]
elif self.kind == "PointInPlane":
markers = [(self._getPnt(args[0]), self._getPln(args[1]))]
elif self.kind == "PointOnLine":
markers = [(self._getPnt(args[0]), self._getLin(args[1]))]
elif self.kind == "Fixed":
markers = [(None,)]
elif self.kind == "FixedPoint":
markers = [(self._getPnt(args[0]),)]
elif self.kind == "FixedAxis":
markers = [(self._getAxis(args[0]),)]
elif self.kind == "FixedRotation":
markers = [(None,), (None,), (None,)]
elif self.kind == "FixedRotationAxis":
markers = [(None,)]
else:
raise ValueError(f"Unknown constraint kind {self.kind}")
# specify kinds of the simple constraint
if self.kind in CompoundConstraints:
kinds, converter = CompoundConstraints[self.kind]
params = converter(self.param,)
else:
kinds = (self.kind,)
params = (self.param,)
# builds the tuple and return
return tuple(zip(markers, kinds, params))
# Cost functions of simple constraints
def point_cost(
m1: gp_Pnt, m2: gp_Pnt, t1: gp_Trsf, t2: gp_Trsf, val: Optional[float] = None,
) -> float:
val = 0 if val is None else val
return val - (m1.Transformed(t1).XYZ() - m2.Transformed(t2).XYZ()).Modulus()
def axis_cost(
m1: gp_Dir, m2: gp_Dir, t1: gp_Trsf, t2: gp_Trsf, val: Optional[float] = None,
) -> float:
val = pi if val is None else val
return DIR_SCALING * (val - m1.Transformed(t1).Angle(m2.Transformed(t2)))
def point_in_plane_cost(
m1: gp_Pnt, m2: gp_Pln, t1: gp_Trsf, t2: gp_Trsf, val: Optional[float] = None,
) -> float:
val = 0 if val is None else val
m2_located = m2.Transformed(t2)
# offset in the plane's normal direction by val:
m2_located.Translate(gp_Vec(m2_located.Axis().Direction()).Multiplied(val))
return m2_located.Distance(m1.Transformed(t1))
def point_on_line_cost(
m1: gp_Pnt, m2: gp_Lin, t1: gp_Trsf, t2: gp_Trsf, val: Optional[float] = None,
) -> float:
val = 0 if val is None else val
m2_located = m2.Transformed(t2)
return val - m2_located.Distance(m1.Transformed(t1))
def fixed_cost(m1: Type[None], t1: gp_Trsf, val: Optional[Type[None]] = None):
return 0
def fixed_point_cost(m1: gp_Pnt, t1: gp_Trsf, val: Tuple[float, float, float]):
return (m1.Transformed(t1).XYZ() - gp_XYZ(*val)).Modulus()
def fixed_axis_cost(m1: gp_Dir, t1: gp_Trsf, val: Tuple[float, float, float]):
return DIR_SCALING * (m1.Transformed(t1).Angle(gp_Dir(*val)))
def fixed_rotation_axis_cost(m1: gp_Dir, t1: gp_Trsf, val: Tuple[int, float]):
ix, v0 = val
v = t1.GetRotation().GetEulerAngles(gp_Intrinsic_XYZ)[ix]
return v - v0
# dictionary of individual constraint cost functions
costs: Dict[str, Callable[..., float]] = dict(
Point=point_cost,
Axis=axis_cost,
PointInPlane=point_in_plane_cost,
PointOnLine=point_on_line_cost,
Fixed=fixed_cost,
FixedPoint=fixed_point_cost,
FixedAxis=fixed_axis_cost,
FixedRotationAxis=fixed_rotation_axis_cost,
)
# Actual solver class
class ConstraintSolver(object):
entities: List[DOF6]
constraints: List[Tuple[Tuple[int, Optional[int]], Constraint]]
constraints: List[Tuple[Tuple[int, ...], Constraint]]
locked: List[int]
ne: int
nc: int
@ -32,24 +402,14 @@ class ConstraintSolver(object):
def __init__(
self,
entities: List[Location],
constraints: List[Tuple[Tuple[int, int], Constraint]],
constraints: List[Tuple[Tuple[int, ...], Constraint]],
locked: List[int] = [],
):
self.entities = [self._locToDOF6(loc) for loc in entities]
self.constraints = []
# decompose into simple constraints
for k, v in constraints:
ms1, ms2, d = v
if ms2:
for m1, m2 in zip(ms1, ms2):
self.constraints.append((k, ((m1,), (m2,), d)))
else:
raise NotImplementedError(
"Single marker constraints are not implemented"
)
self.constraints = constraints
# additional book-keeping
self.ne = len(entities)
self.locked = locked
self.nc = len(self.constraints)
@ -91,82 +451,34 @@ class ConstraintSolver(object):
Callable[[Array[(Any,), float]], float],
Callable[[Array[(Any,), float], Array[(Any,), float]], None],
]:
def pt_cost(
m1: gp_Pnt,
m2: gp_Pnt,
t1: gp_Trsf,
t2: gp_Trsf,
val: Optional[float] = None,
) -> float:
val = 0 if val is None else val
return val - (m1.Transformed(t1).XYZ() - m2.Transformed(t2).XYZ()).Modulus()
def dir_cost(
m1: gp_Dir,
m2: gp_Dir,
t1: gp_Trsf,
t2: gp_Trsf,
val: Optional[float] = None,
) -> float:
val = pi if val is None else val
return DIR_SCALING * (val - m1.Transformed(t1).Angle(m2.Transformed(t2)))
def pnt_pln_cost(
m1: gp_Pnt,
m2: gp_Pln,
t1: gp_Trsf,
t2: gp_Trsf,
val: Optional[float] = None,
) -> float:
val = 0 if val is None else val
m2_located = m2.Transformed(t2)
# offset in the plane's normal direction by val:
m2_located.Translate(gp_Vec(m2_located.Axis().Direction()).Multiplied(val))
return m2_located.Distance(m1.Transformed(t1))
constraints = self.constraints
ne = self.ne
delta = DIFF_EPS * eye(NDOF)
def f(x):
"""
Function to be minimized
"""
constraints = self.constraints
ne = self.ne
rv = 0
transforms = [
self._build_transform(*x[NDOF * i : NDOF * (i + 1)]) for i in range(ne)
]
for i, ((k1, k2), (ms1, ms2, d)) in enumerate(constraints):
t1 = transforms[k1] if k1 not in self.locked else gp_Trsf()
t2 = transforms[k2] if k2 not in self.locked else gp_Trsf()
for ks, (ms, kind, params) in constraints:
ts = tuple(
transforms[k] if k not in self.locked else gp_Trsf() for k in ks
)
cost = costs[kind]
for m1, m2 in zip(ms1, ms2):
if isinstance(m1, gp_Pnt) and isinstance(m2, gp_Pnt):
rv += pt_cost(m1, m2, t1, t2, d) ** 2
elif isinstance(m1, gp_Dir):
rv += dir_cost(m1, m2, t1, t2, d) ** 2
elif isinstance(m1, gp_Pnt) and isinstance(m2, gp_Pln):
rv += pnt_pln_cost(m1, m2, t1, t2, d) ** 2
else:
raise NotImplementedError(f"{m1,m2}")
rv += cost(*ms, *ts, params) ** 2
return rv
def grad(x, rv):
constraints = self.constraints
ne = self.ne
delta = DIFF_EPS * eye(NDOF)
rv[:] = 0
transforms = [
@ -179,60 +491,25 @@ class ConstraintSolver(object):
for j in range(NDOF)
]
for i, ((k1, k2), (ms1, ms2, d)) in enumerate(constraints):
t1 = transforms[k1] if k1 not in self.locked else gp_Trsf()
t2 = transforms[k2] if k2 not in self.locked else gp_Trsf()
for ks, (ms, kind, params) in constraints:
ts = tuple(
transforms[k] if k not in self.locked else gp_Trsf() for k in ks
)
cost = costs[kind]
for m1, m2 in zip(ms1, ms2):
if isinstance(m1, gp_Pnt) and isinstance(m2, gp_Pnt):
tmp = pt_cost(m1, m2, t1, t2, d)
tmp_0 = cost(*ms, *ts, params)
for j in range(NDOF):
for ix, k in enumerate(ks):
if k in self.locked:
continue
t1j = transforms_delta[k1 * NDOF + j]
t2j = transforms_delta[k2 * NDOF + j]
for j in range(NDOF):
tkj = transforms_delta[k * NDOF + j]
if k1 not in self.locked:
tmp1 = pt_cost(m1, m2, t1j, t2, d)
rv[k1 * NDOF + j] += 2 * tmp * (tmp1 - tmp) / DIFF_EPS
ts_kj = ts[:ix] + (tkj,) + ts[ix + 1 :]
tmp_kj = cost(*ms, *ts_kj, params)
if k2 not in self.locked:
tmp2 = pt_cost(m1, m2, t1, t2j, d)
rv[k2 * NDOF + j] += 2 * tmp * (tmp2 - tmp) / DIFF_EPS
elif isinstance(m1, gp_Dir):
tmp = dir_cost(m1, m2, t1, t2, d)
for j in range(NDOF):
t1j = transforms_delta[k1 * NDOF + j]
t2j = transforms_delta[k2 * NDOF + j]
if k1 not in self.locked:
tmp1 = dir_cost(m1, m2, t1j, t2, d)
rv[k1 * NDOF + j] += 2 * tmp * (tmp1 - tmp) / DIFF_EPS
if k2 not in self.locked:
tmp2 = dir_cost(m1, m2, t1, t2j, d)
rv[k2 * NDOF + j] += 2 * tmp * (tmp2 - tmp) / DIFF_EPS
elif isinstance(m1, gp_Pnt) and isinstance(m2, gp_Pln):
tmp = pnt_pln_cost(m1, m2, t1, t2, d)
for j in range(NDOF):
t1j = transforms_delta[k1 * NDOF + j]
t2j = transforms_delta[k2 * NDOF + j]
if k1 not in self.locked:
tmp1 = pnt_pln_cost(m1, m2, t1j, t2, d)
rv[k1 * NDOF + j] += 2 * tmp * (tmp1 - tmp) / DIFF_EPS
if k2 not in self.locked:
tmp2 = pnt_pln_cost(m1, m2, t1, t2j, d)
rv[k2 * NDOF + j] += 2 * tmp * (tmp2 - tmp) / DIFF_EPS
else:
raise NotImplementedError(f"{m1,m2}")
rv[k * NDOF + j] += 2 * tmp_0 * (tmp_kj - tmp_0) / DIFF_EPS
return f, grad

View File

@ -619,7 +619,7 @@ Plane
=====
The Plane constraint is simply a combination of both the Point and Axis constraints. It is a
convienient shortcut for a commonly used combination of constraints. It can be used to shorten the
convenient shortcut for a commonly used combination of constraints. It can be used to shorten the
previous example from the two constraints to just one:
.. code-block:: diff
@ -640,10 +640,9 @@ previous example from the two constraints to just one:
The result of this code is identical to the above two constraint example.
For the cost function of Plane, please see the Point and Axis sections. The ``param`` argument is
copied into both constraints and should be left as the default value of ``None`` for a "mate" style
For the cost function of Plane, please see the Point and Axis sections. The ``param`` argument is applied to Axis and should be left as the default value for a "mate" style
constraint (two surfaces touching) or can be set to ``0`` for a through surface constraint (see
desciption in the Axis constraint section).
description in the Axis constraint section).
PointInPlane
@ -710,6 +709,178 @@ Where:
show_object(assy)
PointOnLine
===========
PointOnLine positions the center of the first object on the line defined by the second object.
The cost function is:
.. math::
( param - \operatorname{dist}(\vec{ c }, l ) )^2
Where:
- :math:`\vec{ c }` is the center of the first argument,
- :math:`l` is a line created from the second object
- :math:`param` is the parameter of the constraint, which defaults to 0,
- :math:`\operatorname{dist}( \vec{ a }, b)` is the distance between point :math:`\vec{ a }` and
line :math:`b`.
.. cadquery::
import cadquery as cq
b1 = cq.Workplane().box(1,1,1)
b2 = cq.Workplane().sphere(0.15)
assy = (
cq.Assembly()
.add(b1,name='b1')
.add(b2, loc=cq.Location(cq.Vector(0,0,4)), name='b2', color=cq.Color('red'))
)
# fix the position of b1
assy.constrain('b1','Fixed')
# b2 on one of the edges of b1
assy.constrain('b2','b1@edges@>>Z and >>Y','PointOnLine')
# b2 on another of the edges of b1
assy.constrain('b2','b1@edges@>>Z and >>X','PointOnLine')
# effectively b2 will be constrained to be on the intersection of the two edges
assy.solve()
show_object(assy)
FixedPoint
==========
FixedPoint fixes the position of the given argument to be equal to the given point specified via the parameter of the constraint. This constraint locks all translational degrees of freedom of the argument.
The cost function is:
.. math::
\left\lVert \vec{ c } - \vec{param} \right\rVert ^2
Where:
- :math:`\vec{ c }` is the center of the argument,
- :math:`param` is the parameter of the constraint - tuple specifying the target position,
- :math:`\operatorname{dist}( \vec{ a }, b)` is the distance between point :math:`\vec{ a }` and
line :math:`b`.
.. cadquery::
import cadquery as cq
b1 = cq.Workplane().box(1,1,1)
b2 = cq.Workplane().sphere(0.15)
assy = (
cq.Assembly()
.add(b1,name='b1')
.add(b2, loc=cq.Location(cq.Vector(0,0,4)), name='b2', color=cq.Color('red'))
)
# fix the position of b1
assy.constrain('b1','Fixed')
# b2 on one of the edges of b1
assy.constrain('b2','b1@edges@>>Z and >>Y','PointOnLine')
# b2 on another of the edges of b1
assy.constrain('b2','b1@edges@>>Z and >>X','PointOnLine')
# effectively b2 will be constrained to be on the intersection of the two edges
assy.solve()
show_object(assy)
FixedRotation
=============
FixedRotation fixes the rotation of the given argument to be equal to the value specified via the parameter of the constraint. The argument is rotated about it's origin firstly by the Z angle, then Y and finally X.
This constraint locks all rotational degrees of freedom of the argument.
The cost function is:
.. math::
\left\lVert \vec{ R } - \vec{param} \right\rVert ^2
Where:
- :math:`\vec{ R }` vector of the rotation angles of the rotation applied to the argument,
- :math:`param` is the parameter of the constraint - tuple specifying the target rotation.
.. cadquery::
import cadquery as cq
b1 = cq.Workplane().box(1,1,1)
b2 = cq.Workplane().rect(0.1, 0.1).extrude(1,taper=-15)
assy = (
cq.Assembly()
.add(b1,name='b1')
.add(b2, loc=cq.Location(cq.Vector(0,0,4)), name='b2', color=cq.Color('red'))
)
# fix the position of b1
assy.constrain('b1','Fixed')
# fix b2 bottom face position (but not rotation)
assy.constrain('b2@faces@<Z','FixedPoint',(0,0,0.5))
# fix b2 rotational degrees of freedom too
assy.constrain('b2','FixedRotation',(45,0,45))
assy.solve()
show_object(assy)
FixedAxis
=========
FixedAxis fixes the orientation of the given argument's normal or tangent to be equal to the orientation of the vector specified via the parameter of the constraint. This constraint locks two rotational degrees of freedom of the argument.
The cost function is:
.. math::
( \vec{ a } \angle \vec{ param } ) ^2
Where:
- :math:`\vec{ a }` normal or tangent vector of the argument,
- :math:`param` is the parameter of the constraint - tuple specifying the target direction.
.. cadquery::
import cadquery as cq
b1 = cq.Workplane().box(1,1,1)
b2 = cq.Workplane().rect(0.1, 0.1).extrude(1,taper=-15)
assy = (
cq.Assembly()
.add(b1,name='b1')
.add(b2, loc=cq.Location(cq.Vector(0,0,4)), name='b2', color=cq.Color('red'))
)
# fix the position of b1
assy.constrain('b1','Fixed')
# fix b2 bottom face position (but not rotation)
assy.constrain('b2@faces@<Z','FixedPoint',(0,0,0.5))
# fix b2 some rotational degrees of freedom too
assy.constrain('b2@faces@>Z','FixedAxis',(1,0,2))
assy.solve()
show_object(assy)
Assembly colors
---------------

View File

@ -106,6 +106,19 @@ def metadata_assy():
return assy
@pytest.fixture
def simple_assy2():
b1 = cq.Workplane().box(1, 1, 1)
b2 = cq.Workplane().box(2, 1, 1)
assy = cq.Assembly()
assy.add(b1, name="b1")
assy.add(b2, loc=cq.Location(cq.Vector(0, 0, 4)), name="b2")
return assy
def test_metadata(metadata_assy):
"""Verify the metadata is present in both the base and sub assemblies"""
assert metadata_assy.metadata["b1"] == "base-data"
@ -327,7 +340,7 @@ def test_constrain(simple_assy, nested_assy):
def test_constrain_with_tags(nested_assy):
nested_assy.add(None, name="dummy")
nested_assy.constrain("TOP?top_face", "SECOND/BOTTOM", "Plane")
nested_assy.constrain("TOP?top_face", "SECOND/BOTTOM", "Point")
assert len(nested_assy.constraints) == 1
@ -445,9 +458,8 @@ def test_constraint_getPln():
return cq.Constraint(ids, (shape0, shape0), sublocs, "PointInPlane", 0)
def fail_this(shape0):
c0 = make_constraint(shape0)
with pytest.raises(ValueError):
c0._getPln(c0.args[0])
make_constraint(shape0)
def resulting_pln(shape0):
c0 = make_constraint(shape0)
@ -593,3 +605,65 @@ def test_infinite_face_constraint_Plane(kind):
)
assy.solve()
assert solve_result_check(assy._solve_result)
def test_unary_constraints(simple_assy2):
assy = simple_assy2
assy.constrain("b1", "Fixed")
assy.constrain("b2", "FixedPoint", (0, 0, -3))
assy.constrain("b2@faces@>Z", "FixedAxis", (0, 1, 1))
assy.solve()
w = cq.Workplane().add(assy.toCompound())
assert w.solids(">Z").val().Center().Length == pytest.approx(0)
assert w.solids("<Z").val().Center().z == pytest.approx(-3)
assert w.solids("<Z").edges(">Z").size() == 1
def test_fixed_rotation(simple_assy2):
assy = simple_assy2
assy.constrain("b1", "Fixed")
assy.constrain("b2", "FixedPoint", (0, 0, -3))
assy.constrain("b2@faces@>Z", "FixedRotation", (45, 0, 0))
assy.solve()
w = cq.Workplane().add(assy.toCompound())
assert w.solids(">Z").val().Center().Length == pytest.approx(0)
assert w.solids("<Z").val().Center().z == pytest.approx(-3)
assert w.solids("<Z").edges(">Z").size() == 1
def test_validation(simple_assy2):
with pytest.raises(ValueError):
simple_assy2.constrain("b1", "Fixed?")
with pytest.raises(ValueError):
cq.assembly.Constraint((), (), (), "Fixed?")
def test_point_on_line(simple_assy2):
assy = simple_assy2
assy.constrain("b1", "Fixed")
assy.constrain("b2@faces@>Z", "FixedAxis", (0, 2, 1))
assy.constrain("b2@faces@>X", "FixedAxis", (1, 0, 0))
assy.constrain("b2@faces@>X", "b1@edges@>>Z and >>Y", "PointOnLine")
assy = assy.solve()
w = cq.Workplane().add(assy.toCompound())
assert w.solids("<Z").val().Center().Length == pytest.approx(0)
assert w.solids(">Z").val().Center().z == pytest.approx(0.5)
assert w.solids(">Z").val().Center().y == pytest.approx(0.5)
assert w.solids(">Z").val().Center().x == pytest.approx(0.0)