Coverage for python/src/dolfinx_mpc/container.py: 95%
79 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-03 12:11 +0000
1from __future__ import annotations
3import warnings
4from typing import Callable, Optional, Union
6import dolfinx
7import numpy
8import numpy.typing as npt
10import dolfinx_mpc.cpp.mpc
12_mpc_data_classes = Union[
13 dolfinx_mpc.cpp.mpc.mpc_data_double,
14 dolfinx_mpc.cpp.mpc.mpc_data_float,
15 dolfinx_mpc.cpp.mpc.mpc_data_complex_double,
16 dolfinx_mpc.cpp.mpc.mpc_data_complex_float,
17]
18_float_array_types = Union[
19 npt.NDArray[numpy.float32],
20 npt.NDArray[numpy.float64],
21 npt.NDArray[numpy.complex64],
22 npt.NDArray[numpy.complex128],
23]
25_mpc_classes = Union[
26 dolfinx_mpc.cpp.mpc.MultiPointConstraint_double,
27 dolfinx_mpc.cpp.mpc.MultiPointConstraint_float,
28 dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_double,
29 dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_float,
30]
31_float_classes = Union[numpy.float32, numpy.float64, numpy.complex128, numpy.complex64]
33_type_names = {
34 numpy.float32: "float",
35 numpy.float64: "double",
36 numpy.complex64: "complex_float",
37 numpy.complex128: "complex_double",
38}
41class _Unset:
42 """The default of a deprecated argument, to tell whether it was passed."""
44 def __repr__(self) -> str:
45 return "<unset>"
48_UNSET = _Unset()
51def _deprecated(old: str, new: str, stacklevel: int = 3):
52 """Warn that the argument `old` is deprecated in favour of `new`, at the caller of the public
53 function `stacklevel - 2` frames up."""
54 warnings.warn(f"`{old}` is deprecated, use {new} instead.", DeprecationWarning, stacklevel=stacklevel)
57def _default_tolerance(dtype: npt.DTypeLike) -> float:
58 """The default distance and coefficient tolerance: 500 machine epsilon of the real type of `dtype`.
60 Mirrors `dolfinx_mpc::default_tolerance` in C++.
61 """
62 return float(500 * numpy.finfo(dtype).eps)
65def _tolerance(value: Optional[float], dtype: npt.DTypeLike) -> float:
66 """`value` as a Python float, by default :func:`_default_tolerance` of `dtype`."""
67 return _default_tolerance(dtype) if value is None else float(value)
70def _cpp_function(name: str, dtype: npt.DTypeLike) -> Callable:
71 """The C++ function `name` for constraints of scalar type `dtype`, bound as `name_<type>`."""
72 return getattr(dolfinx_mpc.cpp.mpc, f"{name}_{_type_names[numpy.dtype(dtype).type]}")
75def _scalar_type(real_type: npt.DTypeLike, dtype: npt.DTypeLike | None = None) -> type:
76 """The scalar type of a constraint on a mesh with coordinates of `real_type`.
78 Args:
79 real_type: The type of the mesh coordinates
80 dtype: The scalar type asked for. Defaults to the default scalar type of DOLFINx, real or
81 complex, at the precision of the mesh.
83 Raises:
84 ValueError: If `dtype` is not one of float32, float64, complex64 and complex128, or its
85 precision differs from the mesh's.
86 """
87 real = numpy.dtype(real_type)
88 if dtype is None:
89 is_complex = numpy.issubdtype(dolfinx.default_scalar_type, numpy.complexfloating)
90 dtype = numpy.promote_types(real, numpy.complex64) if is_complex else real
91 scalar = numpy.dtype(dtype)
92 if scalar.type not in (numpy.float32, numpy.float64, numpy.complex64, numpy.complex128):
93 raise ValueError(f"Unsupported scalar type {scalar} for a constraint")
94 if numpy.finfo(scalar).dtype != real:
95 raise ValueError(
96 f"A constraint of scalar type {scalar} needs a mesh of {numpy.finfo(scalar).dtype}, not {real}"
97 )
98 return scalar.type
101class MPCData:
102 r"""The rows of a constraint on one process, before it is finalized.
104 Row `i` is the equation :math:`u_s = \sum_j c_j u_{m_j}` of the slave :math:`s` =
105 `slaves[i]`, with :math:`m_j` = `masters[j]` and :math:`c_j` = `coeffs[j]` for
106 `offsets[i] <= j < offsets[i + 1]`.
108 Attributes:
109 slaves: The slave of each row: a dof of the space of the slaves, local to the process
110 (owned or ghost) and unrolled
111 masters: The masters of all rows, row after row: global, unrolled dofs of the space of
112 each master
113 coeffs: The coefficient of each master
114 owners: The process owning each master
115 offsets: The masters of row `i` are `masters[offsets[i]:offsets[i + 1]]`, so `offsets`
116 has one entry more than `slaves`
117 master_blocks: The block of each master, in the numbering of whoever made the rows, or
118 `None` if every master is in the space of the slaves. Given on every process or on
119 none. A :class:`MultiPointConstraint` stores codes here until it is finalized (see
120 :meth:`MultiPointConstraint.__init__`).
122 The arrays are views of the C++ object, valid while this object is.
123 """
125 _cpp_object: _mpc_data_classes
127 def __init__(
128 self,
129 slaves: npt.NDArray[numpy.int32],
130 masters: npt.NDArray[numpy.int64],
131 coeffs: _float_array_types,
132 owners: npt.NDArray[numpy.int32],
133 offsets: npt.NDArray[numpy.int32],
134 master_blocks: Optional[npt.NDArray[numpy.int32]] = None,
135 ):
136 args = (
137 numpy.asarray(slaves, dtype=numpy.int32),
138 numpy.asarray(masters, dtype=numpy.int64),
139 coeffs,
140 numpy.asarray(owners, dtype=numpy.int32),
141 numpy.asarray(offsets, dtype=numpy.int32),
142 None if master_blocks is None else numpy.asarray(master_blocks, dtype=numpy.int32),
143 )
144 if coeffs.dtype.type == numpy.float32:
145 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_float(*args)
146 elif coeffs.dtype.type == numpy.float64:
147 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_double(*args)
148 elif coeffs.dtype.type == numpy.complex64:
149 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_complex_float(*args)
150 elif coeffs.dtype.type == numpy.complex128:
151 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_complex_double(*args)
152 else:
153 raise ValueError(f"Unsupported dtype {coeffs.dtype.type} for coefficients")
155 @classmethod
156 def empty(cls, dtype: npt.DTypeLike, master_blocks: bool = False) -> MPCData:
157 """No rows, with coefficients of `dtype`, and with (empty) blocks if `master_blocks`."""
158 return cls(
159 numpy.zeros(0, dtype=numpy.int32),
160 numpy.zeros(0, dtype=numpy.int64),
161 numpy.zeros(0, dtype=dtype),
162 numpy.zeros(0, dtype=numpy.int32),
163 numpy.zeros(1, dtype=numpy.int32),
164 numpy.zeros(0, dtype=numpy.int32) if master_blocks else None,
165 )
167 @classmethod
168 def from_cpp(cls, cpp_object: _mpc_data_classes) -> MPCData:
169 """Wrap rows made in C++, such as the result of a constraint generator."""
170 data = cls.__new__(cls)
171 data._cpp_object = cpp_object
172 return data
174 def append(
175 self,
176 slaves: npt.NDArray[numpy.int32],
177 masters: npt.NDArray[numpy.int64],
178 coeffs: _float_array_types,
179 owners: npt.NDArray[numpy.int32],
180 offsets: npt.NDArray[numpy.int32],
181 master_blocks: Optional[npt.NDArray[numpy.int32]] = None,
182 ) -> MPCData:
183 """These rows followed by the rows given, as new data with the coefficient type of these.
185 The blocks are kept if either has them; missing ones must then be given by the caller.
186 """
187 if (self.master_blocks is None) != (master_blocks is None):
188 raise ValueError("Either both or neither of the rows must have master blocks")
189 dtype = self.coeffs.dtype
190 return MPCData(
191 numpy.concatenate([self.slaves, slaves]),
192 numpy.concatenate([self.masters, masters]),
193 numpy.concatenate([self.coeffs, numpy.asarray(coeffs, dtype=dtype)]),
194 numpy.concatenate([self.owners, owners]),
195 numpy.concatenate([self.offsets, numpy.asarray(offsets[1:]) + self.offsets[-1]]),
196 None if master_blocks is None else numpy.concatenate([self.master_blocks, master_blocks]),
197 )
199 @property
200 def slaves(self):
201 return self._cpp_object.slaves
203 @property
204 def masters(self):
205 return self._cpp_object.masters
207 @property
208 def coeffs(self):
209 return self._cpp_object.coeffs
211 @property
212 def owners(self):
213 return self._cpp_object.owners
215 @property
216 def offsets(self):
217 return self._cpp_object.offsets
219 @property
220 def master_blocks(self):
221 return self._cpp_object.master_blocks