Coverage for python/src/dolfinx_mpc/multipointconstraint.py: 90%
417 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-03 12:11 +0000
1# Copyright (C) 2020-2026 Jørgen S. Dokken
2#
3# This file is part of DOLFINX_MPC
4#
5# SPDX-License-Identifier: MIT
6from __future__ import annotations
8from typing import Callable, Dict, List, Optional, Sequence, Tuple, Union
10from mpi4py import MPI as _MPI
11from petsc4py import PETSc as _PETSc
13import dolfinx.cpp as _cpp
14import dolfinx.fem as _fem
15import dolfinx.mesh as _mesh
16import numpy
17import numpy.typing as npt
18import ufl
20import dolfinx_mpc.cpp
22from .container import (
23 _UNSET,
24 MPCData,
25 _cpp_function,
26 _deprecated,
27 _float_array_types,
28 _float_classes,
29 _mpc_classes,
30 _mpc_data_classes,
31 _scalar_type,
32 _tolerance,
33 _Unset,
34)
35from .dictcondition import create_dictionary_constraint
36from .integralcondition import create_integral_constraint
37from .rbe import create_rbe2, create_rbe3
40class MultiPointConstraint:
41 """
42 Hold data for multi point constraint relation ships,
43 including new index maps for local assembly of matrices and vectors.
45 The constraint is affine, :math:`x = K x_{red} + g`, where :math:`g` is
46 supplied through `rhs_coeffs` and through the Dirichlet conditions in
47 `bcs`. With neither, :math:`g=0` and the constraint is the usual linear
48 one.
50 Args:
51 V: The function space
52 dtype: The scalar type of the coefficients, used by everything built from this
53 constraint. Defaults to the default scalar type of DOLFINx, real or complex, at the
54 precision of the mesh of `V`. Its precision must be that of the mesh.
55 bcs: Dirichlet boundary conditions for the problem. A master degree of
56 freedom that is constrained by one of these is removed from the
57 equation of its slave, and its contribution folded into the
58 constraint offset :math:`g`. As the offset is recomputed from the
59 current values of the conditions by :func:`update_constants`, time
60 dependent boundary data is supported.
61 rhs_coeffs: Function holding an additional inhomogeneity :math:`g_s`
62 for the slave degrees of freedom, i.e.
63 :math:`u_s = \\sum_j c_j u_{m_j} + g_s`.
64 """
66 _data: MPCData
67 _master_spaces: List[_fem.FunctionSpace]
68 _bcs: List[_fem.DirichletBC]
69 _rhs_coeffs: Optional[_fem.Function]
70 _scale_function: Optional[_fem.Function]
71 _rbe2: List[list]
72 _rbe3: List[tuple]
73 _rbe3_data: Optional[tuple]
74 V: _fem.FunctionSpace
75 _input_space: _fem.FunctionSpace
76 finalized: bool
77 _cpp_object: _mpc_classes
78 _dtype: npt.DTypeLike
79 __slots__ = tuple(__annotations__)
81 def __init__(
82 self,
83 V: _fem.FunctionSpace,
84 dtype: npt.DTypeLike | None = None,
85 bcs: Optional[List[_fem.DirichletBC]] = None,
86 rhs_coeffs: Optional[_fem.Function] = None,
87 ):
88 dtype = _scalar_type(V.mesh.geometry.x.dtype, dtype)
89 # The rows on this process until finalize. Their master_blocks hold a code per master, as
90 # the blocks are known only when the constraints are finalized together:
91 # code >= 0 the block itself, its position among the constraints finalized together
92 # code == -1 the block of this constraint: a master in its own space
93 # code == -2 - s the block of the space self._master_spaces[s]
94 # finalize_multipointconstraints resolves the codes to blocks.
95 self._data = MPCData.empty(dtype, master_blocks=True)
96 # The spaces of masters outside the space of this constraint, as given to add_constraint,
97 # add_constraint_from_mpc_data or extend_masters (master_space), each once, in the order
98 # first given. Every process appends them in the same order, also without local slaves,
99 # so that a code means the same space on every process. A space must be the space of one
100 # of the constraints finalized together with this one, or an uncollapsed subspace of it.
101 self._master_spaces = []
102 self._bcs = [] if bcs is None else list(bcs)
103 if rhs_coeffs is not None:
104 if not rhs_coeffs.x.array.dtype == dtype:
105 raise ValueError("rhs_coeffs must have the same dtype as the MPC")
106 if rhs_coeffs.function_space != V:
107 raise ValueError("rhs_coeffs must be a Function in the space of the constraint")
108 self._rhs_coeffs = rhs_coeffs
109 self._scale_function = None
110 # Per space on a spider mesh: [W, the tied space, the block of W], the block set by finalize
111 self._rbe2 = []
112 # Feet of RBE3 constraints on this space, (V, dofs, spiders, weights) per call, built by
113 # finalize, which keeps the arrays for update_rbe3
114 self._rbe3 = []
115 self._rbe3_data = None
116 self.V = V
117 # Kept after finalize replaces `V` by the extended space, which contains no Dirichlet condition
118 self._input_space = V
119 self.finalized = False
120 self._dtype = dtype
122 def add_constraint(
123 self,
124 V: _fem.FunctionSpace,
125 slaves: npt.NDArray[numpy.int32],
126 masters: npt.NDArray[numpy.int64],
127 coeffs: _float_array_types,
128 owners: npt.NDArray[numpy.int32],
129 offsets: npt.NDArray[numpy.int32],
130 master_space: Optional[_fem.FunctionSpace] = None,
131 master_blocks: Optional[npt.NDArray[numpy.int32]] = None,
132 ):
133 """
134 Add new constraint given by numpy arrays.
136 Args:
137 V: The function space for the constraint
138 slaves: List of all slave dofs (using local dof numbering) on this process
139 masters: List of all master dofs (using global dof numbering) on this process
140 coeffs: The coefficients corresponding to each master.
141 owners: The process each master is owned by.
142 offsets: Array indicating the location in the masters array for the i-th slave
143 in the slaves arrays, i.e.
145 .. highlight:: python
146 .. code-block:: python
148 masters_of_owned_slave[i] = masters[offsets[i]:offsets[i+1]]
150 master_space: The function space all masters belong to, if not `V`. It must be the
151 space of another constraint finalized together with this one by
152 :func:`finalize_multipointconstraints`, or an uncollapsed subspace of it, and `masters` is in
153 the global numbering of that constraint's space.
154 The masters of a slave may then be in another block of a blocked problem.
155 master_blocks: The block of each master, for masters from several spaces: its position
156 in the list of constraints given to :func:`finalize_multipointconstraints`. Each
157 master is in the global numbering of its block. Exclusive with `master_space`.
159 Note:
160 Collective when `master_space` or `master_blocks` is given: every process must call
161 it with the same `master_space`, or with `master_blocks` (possibly empty).
162 """
163 assert V == self.V
164 self._raise_if_finalized()
165 if master_space is not None and master_blocks is not None:
166 raise ValueError("Give either master_space or master_blocks, not both")
167 if master_blocks is not None and len(master_blocks) != len(masters):
168 raise ValueError("master_blocks must have one entry per master")
170 if master_blocks is not None:
171 codes = numpy.asarray(master_blocks, dtype=numpy.int32)
172 elif master_space is not None:
173 codes = numpy.full(len(masters), self._space_code(master_space), dtype=numpy.int32)
174 else:
175 codes = numpy.full(len(masters), -1, dtype=numpy.int32)
176 if len(slaves) > 0:
177 self._data = self._data.append(slaves, masters, coeffs, owners, offsets, codes)
179 def _space_code(self, space: _fem.FunctionSpace) -> int:
180 """The code of the block of the masters in `space`, recorded on every process in the same
181 order so that every process resolves it the same way when the constraints are finalized."""
182 for s, other in enumerate(self._master_spaces):
183 if other is space:
184 return -2 - s
185 self._master_spaces.append(space)
186 return -1 - len(self._master_spaces)
188 def extend_masters(
189 self,
190 slaves: npt.NDArray[numpy.int32],
191 masters: npt.NDArray[numpy.int64],
192 coeffs: _float_array_types,
193 owners: npt.NDArray[numpy.int32],
194 master_space: Optional[_fem.FunctionSpace] = None,
195 ):
196 r"""Add a master to the rows of existing slaves.
198 A row is the equation :math:`u_s = \sum_j c_j u_{m_j}` of its slave :math:`u_s`, whichever
199 constraint made it. Entry `i` adds a term :math:`k u_d` to its right-hand side, giving
200 :math:`u_s = \sum_j c_j u_{m_j} + k u_d`, with :math:`s` = `slaves[i]`,
201 :math:`u_d` = `masters[i]` and :math:`k` = `coeffs[i]`. If :math:`u_d` is already a master
202 of the row, in the same space, :math:`k` is added to its coefficient. The ghost copies of
203 the rows are extended too.
205 A call adds one term per slave. Further terms of the same slave are added by further calls,
206 each with its own `master_space`, so that a row may have masters in several spaces.
208 For instance, a periodic condition :math:`u(x) = u(\mathrm{relation}(x))` extended, by one
209 call per entry :math:`\bar H_{bk}` of a macroscopic gradient on a point mesh, with
210 coefficient :math:`d_k = (x - \mathrm{relation}(x))_k` for component :math:`b`, ties `u` to
211 :math:`\bar H x` plus a periodic field.
213 Args:
214 slaves: The slave of each entry: owned dofs of the space of the constraint (local,
215 unrolled), each already a slave and each at most once.
216 masters: The master of each entry, in the global (unrolled) numbering of `master_space`
217 coeffs: The coefficient of each entry, of the scalar type of the constraint
218 owners: The process owning each master
219 master_space: The space of the masters if not that of the constraint. It must be the
220 space of another constraint finalized together with this one by
221 :func:`finalize_multipointconstraints`, or an uncollapsed subspace of it.
223 Raises:
224 ValueError: On every process, if a slave is not owned, appears twice or has no row yet,
225 or the arrays differ in length.
227 Note:
228 Collective. Must be called by every process, with the same `master_space`, before the
229 constraint is finalized.
230 """
231 self._raise_if_finalized()
232 code = -1 if master_space is None else self._space_code(master_space)
233 # Raises ValueError (as the C++ throws std::invalid_argument), identically on every process
234 new_data = _cpp_function("extend_mpc_data", self._dtype)(
235 self._data._cpp_object,
236 -1,
237 numpy.asarray(slaves, dtype=numpy.int32),
238 numpy.asarray(masters, dtype=numpy.int64),
239 numpy.asarray(coeffs, dtype=self._dtype),
240 numpy.asarray(owners, dtype=numpy.int32),
241 code,
242 self.V._cpp_object,
243 )
244 self._data = MPCData.from_cpp(new_data)
246 def add_master_from_point(
247 self,
248 V: _fem.FunctionSpace,
249 slave_marker: Callable[[numpy.ndarray], numpy.ndarray],
250 coefficient: Union[Callable[[numpy.ndarray], numpy.ndarray], float, complex],
251 point: npt.ArrayLike,
252 distance_tol: Optional[float] = None,
253 ):
254 r"""Add the dofs at a point as a master of the rows of marked slaves.
256 Each slave of `V` at a coordinate :math:`x` marked by `slave_marker` gets the term
257 :math:`k(x)\, u(p)` added to the right-hand side of its row, by :meth:`extend_masters`,
258 with :math:`p` = `point` and :math:`k` = `coefficient`. For a blocked `V`, component
259 :math:`b` of a slave gets component :math:`b` of the dofs at :math:`p`; for a single
260 component pass the uncollapsed subspace, for instance `V.sub(1)`.
262 For instance, the periodic cell of the homogenization demos ties the right edge to the left
263 edge plus the corner :math:`B`, :math:`u(L, y) = u(0, y) + u^B`: the periodic constraint
264 gives :math:`u(L, y) = u(0, y)`, and this function, with the right edge marked,
265 coefficient `1` and `point` :math:`B`, adds :math:`u^B`.
267 Args:
268 V: The space of this constraint, or an uncollapsed subspace of it
269 slave_marker: Marks the coordinates of the slaves, `(3, n)` to `(n,)`. Each marked dof
270 must already be a slave, and the dofs at `point` must not be marked.
271 coefficient: The coefficient :math:`k` of each slave, from its coordinates, `(3, n)`
272 to `(n,)`, or one value for all.
273 point: The point of the master, of `gdim` or 3 coordinates. Exactly one block of dofs
274 of `V` must be there (see :func:`dofs_at_point`).
275 distance_tol: The largest distance between `point` and its dofs. Defaults to `500`
276 machine epsilon of the coordinate type of the mesh.
278 Note:
279 Collective. Must be called by every process, before the constraint is finalized.
280 """
281 dofs_marked, components, x_slaves, _ = _marked_dofs(V, slave_marker)
282 owned = dofs_marked < self.V.dofmap.index_map.size_local * self.V.dofmap.index_map_bs
283 slaves, components, x_slaves = dofs_marked[owned], components[owned], x_slaves[owned]
285 if callable(coefficient):
286 coeffs = numpy.asarray(coefficient(x_slaves.T.copy()), dtype=self._dtype).reshape(-1)
287 else:
288 coeffs = numpy.full(len(slaves), coefficient, dtype=self._dtype)
289 dofs, owner = dofs_at_point(V, point, distance_tol)
290 self.extend_masters(
291 slaves.astype(numpy.int32),
292 dofs[components].astype(numpy.int64),
293 coeffs,
294 numpy.full(len(slaves), owner, dtype=numpy.int32),
295 )
297 def add_integral_constraint(
298 self,
299 weight_form,
300 value,
301 bcs: Optional[List[_fem.DirichletBC]] = None,
302 rtol: numpy.floating | float | None = None,
303 *,
304 coefficient_tol: Optional[float] = None,
305 ):
306 r"""Constrain a scalar integral of the solution, :math:`L(u) = \gamma`.
308 The functional is given as a linear form, and turned into a constraint
309 with a single slave by :func:`dolfinx_mpc.create_integral_constraint`;
310 see there for the derivation and the cost. The inhomogeneity
311 :math:`\gamma/w_s` is written into the ``rhs_coeffs`` function of this
312 constraint, which is created here if none was supplied to the
313 constructor.
315 Args:
316 weight_form: A linear form in ``ufl.TestFunction(V)`` defining the
317 functional, for instance ``v * ufl.dx``. Its test function must
318 be in the function space of this constraint.
319 value: The prescribed value :math:`\gamma` of the functional.
320 bcs: Dirichlet conditions on the space. A constrained degree of
321 freedom is never chosen as the slave. Pass the same conditions
322 to the constructor to have a constrained *master* folded into
323 the constraint offset. Defaults to the conditions given to the
324 constructor.
325 rtol: Deprecated, use `coefficient_tol`.
326 coefficient_tol: A master whose coefficient is below
327 `coefficient_tol` times the largest one is dropped. Defaults to
328 `500` machine epsilon of the real type of the constraint.
330 Note:
331 Collective. Must be called by every process.
332 """
333 self._raise_if_finalized()
334 if rtol is not None:
335 _deprecated("rtol", "`coefficient_tol`")
336 coefficient_tol = rtol if coefficient_tol is None else coefficient_tol
337 slaves, masters, coeffs, owners, offsets, rhs = create_integral_constraint(
338 self.V,
339 weight_form,
340 value,
341 self._bcs if bcs is None else bcs,
342 coefficient_tol=_tolerance(coefficient_tol, self._dtype),
343 )
344 if self._rhs_coeffs is None:
345 self._rhs_coeffs = rhs
346 else:
347 # Slaves of separate constraints are disjoint, so the offsets add
348 self._rhs_coeffs.x.array[:] += rhs.x.array
349 self.add_constraint(self.V, slaves, masters, coeffs, owners, offsets)
351 def add_constraint_from_mpc_data(
352 self,
353 V: _fem.FunctionSpace,
354 mpc_data: Union[_mpc_data_classes, MPCData],
355 master_space: Optional[_fem.FunctionSpace] = None,
356 ):
357 """
358 Add new constraint given by an `dolfinc_mpc.cpp.mpc.mpc_data`-object. See
359 :meth:`add_constraint` for `master_space`. The `master_blocks` of `mpc_data`, if any,
360 are the `master_blocks` of :meth:`add_constraint`, and exclusive with `master_space`.
361 """
362 self._raise_if_finalized()
363 self.add_constraint(
364 V,
365 mpc_data.slaves,
366 mpc_data.masters,
367 mpc_data.coeffs,
368 mpc_data.owners,
369 mpc_data.offsets,
370 master_space=master_space,
371 master_blocks=mpc_data.master_blocks,
372 )
374 def finalize(self, filter: Optional[numpy.floating] = None) -> None:
375 """
376 Finializes the multi point constraint. After this function is called, no new constraints can be added
377 to the constraint. This function creates a map from the cells (local to index) to the slave degrees of
378 freedom and builds a new index map and function space where unghosted master dofs are added as ghosts.
380 Args:
381 filter: If given, discard every master whose coefficient satisfies
382 :math:`|c_{sj}| < \\mathrm{filter}\\cdot\\max_k|c_{sk}|`, the
383 maximum being over the masters of that same slave. A negligible
384 coefficient contributes nothing to the constraint, but still
385 costs a ghost, a row of the sparsity pattern and an entry in
386 every element matrix modification, so removing them can shrink
387 :math:`K^HAK` substantially. With `None` (the default) every
388 master supplied is kept.
390 Note:
391 Filtering changes the constraint that is enforced, by exactly the
392 terms that are dropped. It is local and adds no communication.
394 Note:
395 To finalize the constraints of several function spaces, for instance the blocks of a
396 :class:`ufl.MixedFunctionSpace`, use :func:`finalize_multipointconstraints`.
397 """
398 finalize_multipointconstraints([self], filter)
400 def update_constants(self) -> None:
401 """
402 Recompute the constraint offset :math:`g` from the current values of the Dirichlet
403 conditions supplied to the constructor.
405 Call this whenever the value of one of those conditions changes, for instance between
406 time steps, before re-assembling. :class:`LinearProblem` calls it automatically.
408 Note:
409 Collective. Must be called by every process.
410 """
411 self._raise_if_not_finalized()
412 if self._rhs_coeffs is not None:
413 # Pass the array natively. Zero-copy, zero-allocation.
414 num_dofs_local = self.V.dofmap.index_map_bs * (
415 self.V.dofmap.index_map.size_local + self.V.dofmap.index_map.num_ghosts
416 )
417 rhs_coeffs = self._rhs_coeffs.x.array[:num_dofs_local]
418 self._cpp_object.set_rhs_coeffs(rhs_coeffs)
420 self._cpp_object.update_constants()
422 @property
423 def constants(self) -> _float_array_types:
424 """
425 The constraint offset :math:`g` for each degree of freedom local to the process,
426 i.e. the affine term in :math:`x = K x_{red} + g`.
427 """
428 self._raise_if_not_finalized()
429 return self._cpp_object.constants
431 @property
432 def has_inhomogeneity(self) -> bool:
433 """
434 Whether any process carries a non-zero constraint offset. The value is globally
435 reduced, so it is identical on every process.
436 """
437 self._raise_if_not_finalized()
438 return self._cpp_object.has_inhomogeneity
440 @property
441 def master_blocks(self) -> npt.NDArray[numpy.int32]:
442 """
443 The block of each master, parallel to ``masters.array``: the position, in the list given
444 to :func:`finalize_multipointconstraints`, of the constraint whose space the master is in.
445 The local index of a master is in the space of its block.
446 """
447 self._raise_if_not_finalized()
448 return self._cpp_object.master_blocks
450 @property
451 def dtype(self) -> type:
452 """The scalar type of the coefficients, also that of everything built from the constraint."""
453 return self._dtype
455 @property
456 def has_cross_block_masters(self) -> bool:
457 """Whether a master on any process is in another block than the slaves."""
458 self._raise_if_not_finalized()
459 return self._cpp_object.has_cross_block_masters
461 def create_periodic_constraint_topological(
462 self,
463 V: _fem.FunctionSpace,
464 meshtag: _mesh.MeshTags,
465 tag: int,
466 relation: Callable[[numpy.ndarray], numpy.ndarray],
467 bcs: List[_fem.DirichletBC],
468 scale: Union[_float_classes, float, complex] = 1.0,
469 tol: Union[_float_classes, float, None, _Unset] = _UNSET,
470 num_threads: Optional[int] = 1,
471 *,
472 distance_tol: Optional[float] = None,
473 coefficient_tol: Optional[float] = None,
474 ):
475 """
476 Create periodic condition for all closure dofs of on all entities in `meshtag` with value `tag`.
477 :math:`u(x_i) = scale * u(relation(x_i))` for all of :math:`x_i` on marked entities.
479 Args:
480 V: The function space to assign the condition to. Should either be the space of the MPC or an uncollapsed
481 subspace of it.
482 meshtag: MeshTag for entity to apply the periodic condition on
483 tag: Tag indicating which entities should be slaves
484 relation: Lambda-function describing the geometrical relation
485 bcs: Dirichlet boundary conditions for the problem (Periodic constraints will be ignored for these dofs)
486 scale: Factor of the masters, of the scalar type of the constraint
487 tol: Deprecated, use `distance_tol` and `coefficient_tol`: a value sets both, `None` sets
488 `coefficient_tol=0`.
489 num_threads: The number of threads to use for certain operations
490 distance_tol: The largest distance from a mapped slave point to a master cell for the point to
491 be in the cell, and the padding of the bounding boxes of the cells. Defaults to `500`
492 machine epsilon of the coordinate type of the mesh.
493 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
494 its slave is dropped. `0` keeps every master, so that the coefficients can later be changed
495 with :func:`scale_coefficients` or :func:`update_coefficients`. Defaults to `500`
496 machine epsilon of the real type of the constraint.
497 """
498 bcs_ = [bc._cpp_object for bc in bcs]
499 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
500 scale = scale.item() # type: ignore
501 distance_tol, coefficient_tol = self._tolerances(distance_tol, coefficient_tol, tol=tol)
502 is_input_space = V is self.V
503 if not (is_input_space or self.V.contains(V)):
504 raise RuntimeError("The input space has to be an uncollapsed subspace (or the full space) of the MPC")
505 mpc_data = _cpp_function("create_periodic_constraint_topological", self._dtype)(
506 V._cpp_object,
507 meshtag._cpp_object,
508 tag,
509 relation,
510 bcs_,
511 scale,
512 not is_input_space,
513 distance_tol,
514 coefficient_tol,
515 num_threads,
516 )
517 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
519 def create_periodic_constraint_geometrical(
520 self,
521 V: _fem.FunctionSpace,
522 indicator: Callable[[numpy.ndarray], numpy.ndarray],
523 relation: Callable[[numpy.ndarray], numpy.ndarray],
524 bcs: List[_fem.DirichletBC],
525 scale: Union[_float_classes, float, complex] = 1.0,
526 tol: Union[_float_classes, float, None, _Unset] = _UNSET,
527 num_threads: Optional[int] = 1,
528 *,
529 distance_tol: Optional[float] = None,
530 coefficient_tol: Optional[float] = None,
531 ):
532 """
533 Create a periodic condition for all degrees of freedom whose physical location satisfies
534 :math:`indicator(x_i)==True`, i.e.
535 :math:`u(x_i) = scale * u(relation(x_i))` for all :math:`x_i`
537 Args:
538 V: The function space to assign the condition to. Should either be the space of the MPC or an uncollapsed
539 subspace of it.
540 indicator: Lambda-function to locate degrees of freedom that should be slaves
541 relation: Lambda-function describing the geometrical relation to master dofs
542 bcs: Dirichlet boundary conditions for the problem
543 (Periodic constraints will be ignored for these dofs)
544 scale: Factor of the masters, of the scalar type of the constraint
545 tol: Deprecated, use `distance_tol` and `coefficient_tol`: a value sets both, `None` sets
546 `coefficient_tol=0`.
547 num_threads: The number of threads to use for certain operations.
548 distance_tol: The largest distance from a mapped slave point to a master cell for the point to
549 be in the cell, and the padding of the bounding boxes of the cells. Defaults to `500`
550 machine epsilon of the coordinate type of the mesh.
551 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
552 its slave is dropped. `0` keeps every master, so that the coefficients can later be changed
553 with :func:`scale_coefficients` or :func:`update_coefficients`. Defaults to `500`
554 machine epsilon of the real type of the constraint.
555 """
556 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
557 scale = scale.item() # type: ignore
558 distance_tol, coefficient_tol = self._tolerances(distance_tol, coefficient_tol, tol=tol)
559 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
560 is_input_space = V is self.V
561 if not (is_input_space or self.V.contains(V)):
562 raise RuntimeError("The input space has to be an uncollapsed subspace (or the full space) of the MPC")
563 mpc_data = _cpp_function("create_periodic_constraint_geometrical", self._dtype)(
564 V._cpp_object,
565 indicator,
566 relation,
567 bcs,
568 scale,
569 not is_input_space,
570 distance_tol,
571 coefficient_tol,
572 num_threads,
573 )
574 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
576 def create_submesh_constraint(
577 self,
578 V: _fem.FunctionSpace,
579 master_space: _fem.FunctionSpace,
580 entity_map: _mesh.EntityMap,
581 bcs: Optional[List[_fem.DirichletBC]] = None,
582 scale: Union[_float_classes, float, complex] = 1.0,
583 tol: Union[_float_classes, float, None, _Unset] = _UNSET,
584 num_threads: int = 1,
585 *,
586 coefficient_tol: Optional[float] = None,
587 ):
588 r"""
589 Tie the degrees of freedom of `V` to `master_space` on a related mesh: a submesh and its
590 parent, related by `entity_map` as returned by :func:`dolfinx.mesh.create_submesh`.
592 Every degree of freedom of `V` in the closure of a cell related to a cell of
593 `master_space` becomes a slave, :math:`u(x_i) = \mathrm{scale}\, u_m(x_i)`, with
594 :math:`u_m` evaluated in the related cell and component `b` tied to component `b`. With `V`
595 on the submesh this is every degree of freedom of `V`, for instance the trace
596 :math:`\bar u = u|_\Gamma` on a submesh of facets; with `V` on the parent it is the
597 degrees of freedom on the submesh. No search is involved: the related cell is a table
598 lookup. For a submesh of facets the parent cell is one attached to the facet, so for a
599 discontinuous `master_space` the side is arbitrary.
601 Args:
602 V: The space of the constraint, or an uncollapsed subspace of it
603 master_space: The space of another constraint finalized together with this one by
604 :func:`finalize_multipointconstraints`, or an uncollapsed subspace of it. Its mesh and the mesh
605 of `V` are the two meshes of `entity_map`, either way round.
606 entity_map: Relates the cells of the submesh to entities of the parent, of
607 codimension 0 or 1
608 bcs: Dirichlet conditions on the space of the constraint. Their degrees of freedom
609 are not made slaves.
610 scale: Factor of the masters, of the scalar type of the constraint
611 tol: Deprecated, use `coefficient_tol`: a value sets it, `None` sets `coefficient_tol=0`.
612 num_threads: The number of threads to use
613 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
614 its slave is dropped. `0` keeps every basis function of the related cell, so that the
615 coefficients can later be changed with :func:`scale_coefficients` or
616 :func:`update_coefficients`. Defaults to `500` machine epsilon of the real type of the
617 constraint. No distance tolerance is needed, as the related cell is not searched for.
619 Raises:
620 ValueError: If `entity_map` does not relate the two meshes, relates entities other
621 than the cells of the submesh, is of codimension above 1, or the spaces have
622 different numbers of components. Raised on every process.
624 Note:
625 Collective.
626 """
627 self._raise_if_finalized()
628 if not (V is self.V or self.V.contains(V)):
629 raise ValueError("V must be the space of the constraint or an uncollapsed subspace of it")
630 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
631 scale = scale.item() # type: ignore
632 if not isinstance(tol, _Unset):
633 _deprecated("tol", "`coefficient_tol`")
634 if coefficient_tol is None:
635 coefficient_tol = 0.0 if tol is None else tol
636 bcs_ = [] if bcs is None else [bc._cpp_object for bc in bcs]
637 mpc_data = _cpp_function("create_submesh_constraint", self._dtype)(
638 V._cpp_object,
639 master_space._cpp_object,
640 entity_map._cpp_object,
641 bcs_,
642 scale,
643 _tolerance(coefficient_tol, self._dtype),
644 num_threads,
645 )
646 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data, master_space=master_space)
648 def _add_rbe2(self, dofs: list[npt.NDArray[numpy.int32]], W: _fem.FunctionSpace, x=None):
649 """Tie `dofs[k]` to spider `k` of `W`, and record `W` for :meth:`update_rbe2`."""
650 spiders = [numpy.full(len(d), k, dtype=numpy.int64) for k, d in enumerate(dofs)]
651 mpc_data = create_rbe2(
652 self.V,
653 numpy.concatenate(dofs) if dofs else numpy.zeros(0, dtype=numpy.int32),
654 numpy.concatenate(spiders) if spiders else numpy.zeros(0, dtype=numpy.int64),
655 W,
656 self._dtype,
657 x,
658 )
659 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data, master_space=W)
660 if not any(W is entry[0] for entry in self._rbe2):
661 self._rbe2.append([W, self.V, None])
663 def add_rbe2_topological(
664 self,
665 dim: int,
666 entities: Union[npt.NDArray[numpy.int32], Sequence[Optional[npt.NDArray[numpy.int32]]]],
667 W: _fem.FunctionSpace,
668 ):
669 r"""
670 Tie the dofs on mesh entities rigidly to a point, as the RBE2 element of other codes
671 (a rigid "spider").
673 Each dof of this constraint's space on the entities is a "foot" of a spider whose "body"
674 is a point of the point mesh of `W`. Every component of a foot follows the motion of its
675 body,
677 .. math::
679 u(x) = t + \theta \times (x - x_c),
681 where :math:`x_c` is the coordinate of the dofs of `W` at the point, :math:`t` its
682 translation and :math:`\theta` its rotation, the dofs of `W` at the point. Without
683 rotations, :math:`u(x) = t`. All rotation terms are kept, also where their coefficient is
684 zero, so :meth:`update_rbe2` can follow the motion of the meshes.
686 Args:
687 dim: Topological dimension of the entities
688 entities: Entities (local to the process) whose dofs are tied to spider 0, or a
689 sequence whose entry `k` holds the entities tied to the spider with input index
690 `k` (see :func:`dolfinx_mpc.create_spider_mesh`). An entry may be `None`.
691 W: Space on the spider mesh (:func:`dolfinx_mpc.create_spider_mesh`). Its value size
692 is the geometric dimension, for translations only, or 6 in 3D and 3 in 2D, for
693 translations and rotations. `W` must be the space of another constraint
694 finalized together with this one by :func:`finalize_multipointconstraints`.
696 Note:
697 Collective. Must be called by every process, with the same number of entries in
698 `entities`.
699 """
700 self._raise_if_finalized()
701 per_spider = [entities] if isinstance(entities, numpy.ndarray) else list(entities)
702 dofs = []
703 # The dofs of each spider in turn, collectively, so that a dof on the entities of two
704 # spiders is caught as constrained twice
705 for e in per_spider:
706 e = numpy.zeros(0, dtype=numpy.int32) if e is None else numpy.asarray(e, dtype=numpy.int32)
707 dofs.append(_fem.locate_dofs_topological(self.V, dim, e))
708 self._add_rbe2(dofs, W)
710 def add_rbe2_geometrical(
711 self,
712 locators: Union[
713 Callable[[numpy.ndarray], numpy.ndarray], Sequence[Optional[Callable[[numpy.ndarray], numpy.ndarray]]]
714 ],
715 W: _fem.FunctionSpace,
716 ):
717 r"""
718 Tie the dofs located geometrically rigidly to a point, as the RBE2 element of other codes
719 (a rigid "spider"). See :meth:`add_rbe2_topological` for the relation.
721 Args:
722 locators: Marks the dofs tied to spider 0, given their coordinates, shape
723 `(3, num_points)`, or a sequence whose entry `k` marks the dofs tied to the spider
724 with input index `k`. An entry may be `None`.
725 W: Space on the spider mesh, see :meth:`add_rbe2_topological`
727 Note:
728 Collective. Must be called by every process, with the same number of locators.
729 """
730 self._raise_if_finalized()
731 per_spider = [locators] if callable(locators) else list(locators)
732 dofs = [
733 numpy.zeros(0, dtype=numpy.int32)
734 if locator is None
735 else numpy.asarray(_fem.locate_dofs_geometrical(self.V, locator), dtype=numpy.int32)
736 for locator in per_spider
737 ]
738 self._add_rbe2(dofs, W, self.V.tabulate_dof_coordinates())
740 def update_rbe2(self) -> None:
741 """
742 Recompute the coefficients of every RBE2 constraint from the current coordinates.
744 The feet are at the dof coordinates of the constraint's space, the spiders at those of the
745 space on the spider mesh, both read now. Move the meshes, for instance to the deformed
746 configuration in an updated Lagrangian analysis, then call this to tie the feet to the
747 rigid motion about the new positions. Assemble again afterwards.
749 The constraint must be finalized without a `filter`, which could drop a master whose
750 coefficient becomes nonzero.
752 Note:
753 Collective. Must be called by every process.
754 """
755 self._raise_if_not_finalized()
756 if len(self._rbe2) == 0:
757 raise ValueError("The constraint has no RBE2 constraints")
758 for W, V, block in self._rbe2:
759 dolfinx_mpc.cpp.mpc.update_rbe2(self._cpp_object, V._cpp_object, W._cpp_object, block)
761 def _add_rbe3(self, V: _fem.FunctionSpace, dofs: list[npt.NDArray[numpy.int32]], weights, x):
762 """Record `dofs[k]` as feet of spider `k`, with weights from `weights`."""
763 real = V.mesh.geometry.x.dtype
764 for k, d in enumerate(dofs):
765 if weights is None:
766 w = numpy.ones(len(d), dtype=real)
767 elif callable(weights):
768 w = numpy.asarray(weights(x[d].T), dtype=real).reshape(-1)
769 else:
770 w = numpy.full(len(d), weights, dtype=real)
771 self._rbe3.append((V, d, numpy.full(len(d), k, dtype=numpy.int64), w))
773 def add_rbe3_topological(
774 self,
775 V: _fem.FunctionSpace,
776 dim: int,
777 entities: Union[npt.NDArray[numpy.int32], Sequence[Optional[npt.NDArray[numpy.int32]]]],
778 weights: Union[None, float, Callable[[numpy.ndarray], numpy.ndarray]] = None,
779 ):
780 r"""
781 Tie the dofs of spiders to the motion of the dofs of `V` on mesh entities, as the RBE3
782 element of other codes (a flexible "spider").
784 This constraint is on the space of the spider mesh (:func:`dolfinx_mpc.create_spider_mesh`).
785 Each spider moves with the rigid motion that best fits its "feet", in the weighted
786 least-squares sense,
788 .. math::
790 \min_{t, \theta} \sum_i w_i |u_i - t - \theta \times (x_i - x_c)|^2,
792 where :math:`x_c` is the coordinate of the spider, :math:`t` and :math:`\theta` its
793 translation and rotation, and :math:`u_i` the displacement of foot :math:`i` at
794 :math:`x_i`. Without rotations, :math:`t` is the weighted mean of the feet. Unlike RBE2, the
795 feet keep their stiffness: a load on the spider is spread over them without making them
796 rigid.
798 The feet may be in several spaces, given by one call each. The constraint is built when
799 it is finalized, by :func:`finalize_multipointconstraints` together with the constraints of
800 the spaces of the feet.
802 Args:
803 V: The space of the feet, with one component per dimension
804 dim: Topological dimension of the entities
805 entities: Entities (local to the process) whose dofs are feet of spider 0, or a
806 sequence whose entry `k` holds those of the spider with input index `k`. An entry
807 may be `None`.
808 weights: The weight of each foot: `None` for one, a number, or a function of the
809 coordinates, shape `(3, num_points)`, returning one non-negative weight per foot.
810 Evaluated once, here: :meth:`update_rbe3` keeps the weights.
812 Note:
813 Collective. Must be called by every process, with the same number of entries in
814 `entities`.
815 """
816 self._raise_if_finalized()
817 per_spider = [entities] if isinstance(entities, numpy.ndarray) else list(entities)
818 dofs = []
819 for e in per_spider:
820 e = numpy.zeros(0, dtype=numpy.int32) if e is None else numpy.asarray(e, dtype=numpy.int32)
821 dofs.append(_fem.locate_dofs_topological(V, dim, e))
822 self._add_rbe3(V, dofs, weights, V.tabulate_dof_coordinates() if callable(weights) else None)
824 def add_rbe3_geometrical(
825 self,
826 V: _fem.FunctionSpace,
827 locators: Union[
828 Callable[[numpy.ndarray], numpy.ndarray], Sequence[Optional[Callable[[numpy.ndarray], numpy.ndarray]]]
829 ],
830 weights: Union[None, float, Callable[[numpy.ndarray], numpy.ndarray]] = None,
831 ):
832 """
833 Tie the dofs of spiders to the motion of the dofs of `V` located geometrically, as the
834 RBE3 element of other codes. See :meth:`add_rbe3_topological` for the relation.
836 Args:
837 V: The space of the feet, with one component per dimension
838 locators: Marks the feet of spider 0, given their coordinates, shape
839 `(3, num_points)`, or a sequence whose entry `k` marks the feet of the spider with
840 input index `k`. An entry may be `None`.
841 weights: See :meth:`add_rbe3_topological`
843 Note:
844 Collective. Must be called by every process, with the same number of locators.
845 """
846 self._raise_if_finalized()
847 per_spider = [locators] if callable(locators) else list(locators)
848 dofs = [
849 numpy.zeros(0, dtype=numpy.int32)
850 if locator is None
851 else numpy.asarray(_fem.locate_dofs_geometrical(V, locator), dtype=numpy.int32)
852 for locator in per_spider
853 ]
854 self._add_rbe3(V, dofs, weights, V.tabulate_dof_coordinates() if callable(weights) else None)
856 def _build_rbe3(self, mpcs: Sequence[MultiPointConstraint]) -> None:
857 """Build the RBE3 constraint from the feet recorded, with the blocks of their spaces."""
858 spaces: list[_fem.FunctionSpace] = []
859 for V, *_ in self._rbe3:
860 if not any(V is other for other in spaces):
861 spaces.append(V)
862 blocks = []
863 for V in spaces:
864 matches = [j for j, other in enumerate(mpcs) if other.V is V]
865 if len(matches) != 1:
866 raise ValueError(
867 "The feet of an RBE3 constraint must be in the space of exactly one of the "
868 "constraints finalized together with it"
869 )
870 blocks.append(matches[0])
872 def gather(i):
873 return [numpy.concatenate([r[i] for r in self._rbe3 if r[0] is V]) for V in spaces]
875 data = (self.V, spaces, gather(1), gather(2), gather(3))
876 mpc_data = create_rbe3(*data, dtype=self._dtype)
877 self.add_constraint(
878 self.V,
879 mpc_data.slaves,
880 mpc_data.masters,
881 mpc_data.coeffs,
882 mpc_data.owners,
883 mpc_data.offsets,
884 master_blocks=numpy.asarray(blocks, dtype=numpy.int32)[mpc_data.master_blocks],
885 )
886 self._rbe3_data = (*data, blocks)
888 def update_rbe3(self) -> None:
889 """
890 Recompute the coefficients of the RBE3 constraint from the current coordinates.
892 The feet are at the dof coordinates of their spaces, the spiders at those of the space of
893 this constraint, both read now. Move the meshes, then call this. Assemble again
894 afterwards.
896 Note:
897 Collective. Must be called by every process.
898 """
899 self._raise_if_not_finalized()
900 if self._rbe3_data is None:
901 raise ValueError("The constraint has no RBE3 constraint")
902 W, spaces, dofs, spiders, weights, blocks = self._rbe3_data
903 dolfinx_mpc.cpp.mpc.update_rbe3(
904 self._cpp_object,
905 W._cpp_object,
906 [V._cpp_object for V in spaces],
907 blocks,
908 [numpy.ascontiguousarray(d, dtype=numpy.int32) for d in dofs],
909 [numpy.ascontiguousarray(k, dtype=numpy.int64) for k in spiders],
910 [numpy.ascontiguousarray(w) for w in weights],
911 )
913 def create_slip_constraint(
914 self,
915 space: _fem.FunctionSpace,
916 facet_marker: Tuple[_mesh.MeshTags, int],
917 v: _fem.Function,
918 bcs: List[_fem.DirichletBC] = [],
919 ):
920 """
921 Create a slip constraint :math:`u \\cdot v=0` over the entities defined in `facet_marker` with the given index.
923 Args:
924 space: Function space (possibly an uncollapsed subspace) for the current constraint
925 facet_marker: Tuple containomg the mesh tag and marker used to locate degrees of freedom
926 v: Function containing the directional vector to dot your slip condition (most commonly a normal vector)
927 bcs: List of Dirichlet BCs (slip conditions will be ignored on these dofs)
929 Examples:
930 Create constaint :math:`u\\cdot n=0` of all indices in `mt` marked with `i`
932 .. highlight:: python
933 .. code-block:: python
935 V = dolfinx.fem.functionspace(mesh, ("CG", 1))
936 mpc = MultiPointConstraint(V)
937 n = dolfinx.fem.Function(V)
938 mpc.create_slip_constaint(V, (mt, i), n)
940 Create slip constaint for a mixed function space:
942 .. highlight:: python
943 .. code-block:: python
945 cellname = mesh.basix_cell()
946 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
947 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
948 me = basix.ufl.mixed_element([Ve, Qe])
949 W = dolfinx.fem.functionspace(mesh, me)
950 mpc = MultiPointConstraint(W)
951 n_space, _ = W.sub(0).collapse()
952 normal = dolfinx.fem.Function(n_space)
953 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[])
955 A slip condition cannot be applied on the same degrees of freedom as a Dirichlet BC, and therefore
956 any Dirichlet bc for the space of the multi point constraint should be supplied.
958 .. highlight:: python
959 .. code-block:: python
961 cellname = mesh.basix_cell()
962 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
963 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
964 me = basix.ufl.mixed_element([Ve, Qe])
965 W = dolfinx.fem.functionspace(mesh, me)
966 mpc = MultiPointConstraint(W)
967 n_space, _ = W.sub(0).collapse()
968 normal = Function(n_space)
969 bc = dolfinx.fem.dirichletbc(inlet_velocity, dofs, W.sub(0))
970 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[bc])
971 """
972 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
973 if space is self.V:
974 sub_space = False
975 elif self.V.contains(space):
976 sub_space = True
977 else:
978 raise ValueError("Input space has to be an uncollapsed subspace of the MPC space")
979 mpc_data = _cpp_function("create_slip_condition", self._dtype)(
980 space._cpp_object,
981 facet_marker[0]._cpp_object,
982 facet_marker[1],
983 v._cpp_object,
984 bcs,
985 sub_space,
986 )
987 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
989 def create_general_constraint(
990 self,
991 slave_master_dict: Dict[bytes, Dict[bytes, float]],
992 subspace_slave: Optional[int] = None,
993 subspace_master: Optional[int] = None,
994 *,
995 distance_tol: Optional[float] = None,
996 ):
997 """
998 Args:
999 V: The function space
1000 slave_master_dict: Nested dictionary, where the first key is the bit representing the slave dof's
1001 coordinate in the mesh. The item of this key is a dictionary, where each key of this dictionary
1002 is the bit representation of the master dof's coordinate, and the item the coefficient for
1003 the MPC equation.
1004 subspace_slave: If using mixed or vector space, and only want to use dofs from a sub space
1005 as slave add index here
1006 subspace_master: Subspace index for mixed or vector spaces
1007 distance_tol: The largest distance between the coordinate of a key and that of its dof.
1008 Defaults to `500` machine epsilon of the coordinate type of the mesh.
1010 Example:
1011 If the dof `D` located at `[d0, d1]` should be constrained to the dofs
1012 `E` and `F` at `[e0, e1]` and `[f0, f1]` as :math:`D = \\alpha E + \\beta F`
1013 the dictionary should be:
1015 .. highlight:: python
1016 .. code-block:: python
1018 {numpy.array([d0, d1], dtype=mesh.geometry.x.dtype).tobytes():
1019 {numpy.array([e0, e1], dtype=mesh.geometry.x.dtype).tobytes(): alpha,
1020 numpy.array([f0, f1], dtype=mesh.geometry.x.dtype).tobytes(): beta}}
1021 """
1022 slaves, masters, coeffs, owners, offsets = create_dictionary_constraint(
1023 self.V,
1024 slave_master_dict,
1025 subspace_slave,
1026 subspace_master,
1027 dtype=self._dtype,
1028 distance_tol=_tolerance(distance_tol, self.V.mesh.geometry.x.dtype),
1029 )
1030 self.add_constraint(self.V, slaves, masters, coeffs, owners, offsets)
1032 def _tolerances(
1033 self,
1034 distance_tol: Optional[float],
1035 coefficient_tol: Optional[float],
1036 tol: Union[_float_classes, float, None, _Unset] = _UNSET,
1037 eps2: Optional[float] = None,
1038 ) -> tuple[float, float]:
1039 """The distance and coefficient tolerance, by default `500` machine epsilon of the coordinate
1040 type of the mesh and of the real type of the constraint.
1042 The deprecated `tol` of the periodic constraints sets both, or with `None` keeps every
1043 master. The deprecated `eps2` of the contact constraints is a squared distance.
1044 """
1045 if not isinstance(tol, _Unset):
1046 _deprecated("tol", "`distance_tol` and `coefficient_tol`", stacklevel=4)
1047 if tol is None:
1048 coefficient_tol = 0.0 if coefficient_tol is None else coefficient_tol
1049 else:
1050 distance_tol = tol if distance_tol is None else distance_tol
1051 coefficient_tol = tol if coefficient_tol is None else coefficient_tol
1052 if eps2 is not None:
1053 _deprecated("eps2", "`distance_tol`, a distance rather than a squared distance,", stacklevel=4)
1054 distance_tol = float(numpy.sqrt(eps2)) if distance_tol is None else distance_tol
1055 return (
1056 _tolerance(distance_tol, self.V.mesh.geometry.x.dtype),
1057 _tolerance(coefficient_tol, self._dtype),
1058 )
1060 def create_contact_slip_condition(
1061 self,
1062 meshtags: _mesh.MeshTags,
1063 slave_marker: int,
1064 master_marker: int,
1065 normal: _fem.Function,
1066 eps2: Optional[float] = None,
1067 num_threads: Optional[int] = 1,
1068 *,
1069 distance_tol: Optional[float] = None,
1070 coefficient_tol: Optional[float] = None,
1071 ):
1072 """
1073 Create a slip condition between two sets of facets marker with individual markers.
1074 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
1075 The condition created is :math:`u_s \\cdot normal_s = u_m \\cdot normal_m` where `s` is the
1076 restriction to the slave facets, `m` to the master facets.
1078 Args:
1079 meshtags: The meshtags of the set of facets to tie together
1080 slave_marker: The marker of the slave facets
1081 master_marker: The marker of the master facets
1082 normal: The function used in the dot-product of the constraint
1083 eps2: Deprecated, use `distance_tol`, which is `sqrt(eps2)`.
1084 num_threads: The number of threads to use for certain operations
1085 distance_tol: The largest distance from a slave point to a master cell for the point to
1086 be in the cell, and the padding of the bounding boxes of the cells. Defaults to `500`
1087 machine epsilon of the coordinate type of the mesh.
1088 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
1089 its slave is dropped. `0` keeps every master. Defaults to `500`
1090 machine epsilon of the real type of the constraint.
1091 """
1092 mpc_data = _cpp_function("create_contact_slip_condition", self._dtype)(
1093 self.V._cpp_object,
1094 meshtags._cpp_object,
1095 slave_marker,
1096 master_marker,
1097 normal._cpp_object,
1098 *self._tolerances(distance_tol, coefficient_tol, eps2=eps2),
1099 num_threads,
1100 )
1101 self.add_constraint_from_mpc_data(self.V, mpc_data)
1103 def create_contact_inelastic_condition(
1104 self,
1105 meshtags: _cpp.mesh.MeshTags_int32,
1106 slave_marker: int,
1107 master_marker: int,
1108 eps2: Optional[float] = None,
1109 allow_missing_masters: bool = False,
1110 num_threads: Optional[int] = 1,
1111 *,
1112 distance_tol: Optional[float] = None,
1113 coefficient_tol: Optional[float] = None,
1114 ):
1115 """
1116 Create a contact inelastic condition between two sets of facets marker with individual markers.
1117 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
1118 The condition created is :math:`u_s = u_m` where `s` is the restriction to the
1119 slave facets, `m` to the master facets.
1121 Args:
1122 meshtags: The meshtags of the set of facets to tie together
1123 slave_marker: The marker of the slave facets
1124 master_marker: The marker of the master facets
1125 eps2: Deprecated, use `distance_tol`, which is `sqrt(eps2)`.
1126 allow_missing_masters: If true, the function will not throw an error if a degree of freedom
1127 in the closure of the master entities does not have a corresponding set of slave degree
1128 of freedom.
1129 num_threads: The number of threads to use for certain operations
1130 distance_tol: The largest distance from a slave point to a master cell for the point to
1131 be in the cell, and the padding of the bounding boxes of the cells. Defaults to `500`
1132 machine epsilon of the coordinate type of the mesh.
1133 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
1134 its slave is dropped. `0` keeps every master. Defaults to `500`
1135 machine epsilon of the real type of the constraint.
1136 """
1137 mpc_data = _cpp_function("create_contact_inelastic_condition", self._dtype)(
1138 self.V._cpp_object,
1139 meshtags._cpp_object,
1140 slave_marker,
1141 master_marker,
1142 *self._tolerances(distance_tol, coefficient_tol, eps2=eps2),
1143 allow_missing_masters,
1144 num_threads,
1145 )
1146 self.add_constraint_from_mpc_data(self.V, mpc_data)
1148 @property
1149 def is_slave(self) -> numpy.ndarray:
1150 """
1151 Returns a vector of integers where the ith entry indicates if a degree of freedom (local to process) is a slave.
1152 """
1153 self._raise_if_not_finalized()
1154 return self._cpp_object.is_slave
1156 @property
1157 def slaves(self):
1158 """
1159 Returns the degrees of freedom for all slaves local to process
1160 """
1161 self._raise_if_not_finalized()
1162 return self._cpp_object.slaves
1164 @property
1165 def masters(self) -> _cpp.graph.AdjacencyList_int32:
1166 """
1167 Returns an adjacency-list whose ith node corresponds to
1168 a degree of freedom (local to process), and links the corresponding master dofs (local to process).
1170 Examples:
1172 .. highlight:: python
1173 .. code-block:: python
1175 masters = mpc.masters
1176 masters_of_dof_i = masters.links(i)
1177 """
1178 self._raise_if_not_finalized()
1179 return self._cpp_object.masters
1181 def coefficients(self) -> _float_array_types:
1182 """
1183 Returns a vector containing the coefficients for the constraint, and the corresponding offsets
1184 for the ith degree of freedom.
1186 Examples:
1188 .. highlight:: python
1189 .. code-block:: python
1191 coeffs, offsets = mpc.coefficients()
1192 coeffs_of_slave_i = coeffs[offsets[i]:offsets[i+1]]
1193 """
1194 self._raise_if_not_finalized()
1195 return self._cpp_object.coefficients()
1197 def all_coefficients(self) -> Tuple[_float_array_types, npt.NDArray[numpy.int32]]:
1198 """
1199 Returns the coefficients of all masters, including those eliminated by a Dirichlet condition,
1200 in the order supplied before :func:`finalize`, and the offsets for the ith degree of freedom.
1201 This is the layout taken by :func:`update_coefficients`. The corresponding masters are given
1202 by :func:`all_masters`.
1204 Examples:
1206 .. highlight:: python
1207 .. code-block:: python
1209 coeffs, offsets = mpc.all_coefficients()
1210 coeffs_of_slave_i = coeffs[offsets[i]:offsets[i+1]]
1211 """
1212 self._raise_if_not_finalized()
1213 return self._cpp_object.all_coefficients()
1215 def all_masters(self) -> npt.NDArray[numpy.int32]:
1216 """
1217 Returns the masters (local index in :attr:`function_space`) in the layout of
1218 :func:`all_coefficients`.
1219 """
1220 self._raise_if_not_finalized()
1221 return self._cpp_object.all_masters()
1223 def update_coefficients(self, coeffs: _float_array_types) -> None:
1224 """
1225 Replace the coefficient of every master, including masters eliminated by a Dirichlet
1226 condition, and recompute the constraint offset :math:`g`.
1228 The masters are fixed at creation. A master dropped by `coefficient_tol` or by the `filter`
1229 of :func:`finalize` cannot be given a coefficient, so create the constraint with
1230 `coefficient_tol=0` and no filter if the coefficients are to be changed.
1232 Args:
1233 coeffs: The new coefficients, in the layout of :func:`all_coefficients`, for all degrees
1234 of freedom local to the process (owned and ghost).
1236 Note:
1237 Collective. Must be called by every process.
1238 """
1239 self._raise_if_not_finalized()
1240 self._cpp_object.update_coefficients(numpy.ascontiguousarray(coeffs, dtype=self._dtype))
1242 def scale_coefficients(
1243 self,
1244 scale: Union[_float_classes, float, complex, ufl.core.expr.Expr, _fem.Expression],
1245 ) -> None:
1246 """
1247 Multiply the coefficients of all masters of each slave :math:`s` by a factor
1248 :math:`f_s`, and recompute the constraint offset :math:`g`. For a periodic constraint
1249 :math:`u(x_s) = f_s u(relation(x_s))`, which for instance gives a Floquet-Bloch condition
1250 with :math:`f=e^{i k\\cdot L}`.
1252 The factors are stored in a function in the space of the constraint, and :math:`f_s` is
1253 the degree of freedom :math:`s` of that function: the value at the slave for a Lagrange
1254 space, the corresponding moment for e.g. a Nédélec space.
1256 Repeated calls compound. Masters eliminated by a Dirichlet condition are scaled as well,
1257 the user supplied `rhs_coeffs` are not.
1259 Args:
1260 scale: A scalar, a :class:`dolfinx.fem.Function` in the constraint's space (copied
1261 by interpolation), a UFL expression, compiled into a :class:`dolfinx.fem.Expression`
1262 at the interpolation points of the space, or such a compiled expression. Pass a
1263 compiled expression to avoid recompilation when the factor is updated through
1264 :class:`dolfinx.fem.Constant`'s in it.
1266 Note:
1267 Collective. Must be called by every process.
1268 """
1269 self._raise_if_not_finalized()
1270 if self._scale_function is None:
1271 self._scale_function = _fem.Function(self.V, dtype=self._dtype)
1272 f = self._scale_function
1273 if isinstance(scale, (_fem.Expression, _fem.Function)):
1274 f.interpolate(scale)
1275 elif isinstance(scale, ufl.core.expr.Expr):
1276 f.interpolate(_fem.Expression(scale, self.V.element.interpolation_points, dtype=self._dtype))
1277 else:
1278 f.x.array[:] = scale
1279 f.x.scatter_forward()
1280 # The extended index map appends master ghosts after the ghosts of the input space
1281 num_dofs_local = len(self._cpp_object.is_slave)
1282 self._cpp_object.scale_coefficients(f.x.array[:num_dofs_local])
1284 @property
1285 def num_local_slaves(self):
1286 """
1287 Return the number of slaves owned by the current process.
1288 """
1289 self._raise_if_not_finalized()
1290 return self._cpp_object.num_local_slaves
1292 @property
1293 def cell_to_slaves(self):
1294 """
1295 Returns an `dolfinx.cpp.graph.AdjacencyList_int32` whose ith node corresponds to
1296 the ith cell (local to process), and links the corresponding slave degrees of
1297 freedom in the cell (local to process).
1299 Examples:
1301 .. highlight:: python
1302 .. code-block:: python
1304 cell_to_slaves = mpc.cell_to_slaves()
1305 slaves_in_cell_i = cell_to_slaves.links(i)
1306 """
1307 self._raise_if_not_finalized()
1308 return self._cpp_object.cell_to_slaves
1310 @property
1311 def function_space(self):
1312 """
1313 Return the function space for the multi-point constraint with the updated index map
1314 """
1315 self._raise_if_not_finalized()
1316 return self.V
1318 @property
1319 def input_space(self) -> _fem.FunctionSpace:
1320 """
1321 The function space the constraint was created with.
1323 Forms and Dirichlet conditions are stated on this space, while functions holding a
1324 solution live in :attr:`function_space`, its extension by the masters of the constraint.
1325 For a system of several blocks, ``[mpc.input_space for mpc in mpcs]`` gives the spaces in
1326 the order of the blocks, for instance for a :class:`ufl.MixedFunctionSpace`.
1327 """
1328 return self._input_space
1330 def backsubstitution(self, u: Union[_fem.Function, Sequence[_fem.Function], _PETSc.Vec]) -> None: # type: ignore
1331 """
1332 For a Function, impose the multi-point constraint by backsubstiution.
1333 This function is used after solving the reduced problem to obtain the values
1334 at the slave degrees of freedom
1336 .. note::
1337 It is the users responsibility to destroy the PETSc vector
1339 Args:
1340 u: The input function. For a constraint with masters in another block, the function
1341 of every block, in the order given to :func:`finalize_multipointconstraints`;
1342 only the function of this constraint's block is changed. The ghosts of the
1343 functions holding masters must be up to date.
1344 """
1345 self._raise_if_not_finalized()
1346 if isinstance(u, Sequence):
1347 self._cpp_object.backsubstitution([u_k.x.array for u_k in u]) # type: ignore
1348 u[self._cpp_object.block].x.scatter_forward()
1349 return
1350 try:
1351 self._cpp_object.backsubstitution(u.x.array) # type: ignore
1352 assert isinstance(u, _fem.Function)
1353 u.x.scatter_forward()
1354 except AttributeError:
1355 assert isinstance(u, _PETSc.Vec)
1356 with u.localForm() as vector_local:
1357 self._cpp_object.backsubstitution(vector_local.array_w)
1358 u.ghostUpdate(addv=_PETSc.InsertMode.INSERT, mode=_PETSc.ScatterMode.FORWARD) # type: ignore
1360 def homogenize(self, u: _fem.Function) -> None:
1361 """
1362 For a vector, homogenize (set to zero) the vector components at the multi-point
1363 constraint slave DoF indices. This is particularly useful for nonlinear problems.
1365 Args:
1366 u: The input vector
1367 """
1368 self._cpp_object.homogenize(u.x.array)
1369 u.x.scatter_forward()
1371 def _raise_if_finalized(self):
1372 """
1373 Raise if the multi point constraint has already been finalized
1374 """
1375 if self.finalized:
1376 raise RuntimeError("MultiPointConstraint has already been finalized")
1378 def _raise_if_not_finalized(self):
1379 """
1380 Raise if the multi point constraint has not yet been finalized
1381 """
1382 if not self.finalized:
1383 raise RuntimeError("MultiPointConstraint has not been finalized")
1386def _marked_dofs(
1387 V: _fem.FunctionSpace, marker: Callable[[numpy.ndarray], numpy.ndarray]
1388) -> tuple[npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.floating], int]:
1389 """The dofs of `V` at the coordinates marked by `marker`, owned and ghosts.
1391 Returns:
1392 Their local, unrolled index in the space `V` is an uncollapsed subspace of (or in `V`),
1393 their component in `V`, their coordinates `(n, 3)`, and the number of components of `V`.
1394 The component is that of `V`, from its collapsed numbering: the block size of the parent
1395 need not be that of `V` (a component of a blocked space, or a block of a mixed space).
1396 """
1397 if len(V.component()) > 0:
1398 V_c = V.collapse()[0]
1399 bs = V_c.dofmap.index_map_bs
1400 parent_dofs, sub_dofs = _fem.locate_dofs_geometrical((V, V_c), marker)
1401 x = V_c.tabulate_dof_coordinates()[sub_dofs // bs]
1402 return parent_dofs, sub_dofs % bs, x, bs
1403 bs = V.dofmap.index_map_bs
1404 blocks = _fem.locate_dofs_geometrical(V, marker)
1405 x = numpy.repeat(V.tabulate_dof_coordinates()[blocks], bs, axis=0)
1406 dofs = (blocks[:, None] * bs + numpy.arange(bs, dtype=numpy.int32)).ravel()
1407 return dofs, numpy.tile(numpy.arange(bs, dtype=numpy.int32), len(blocks)), x, bs
1410def dofs_at_point(
1411 V: _fem.FunctionSpace, point: npt.ArrayLike, distance_tol: float | None = None
1412) -> tuple[npt.NDArray[numpy.int64], int]:
1413 """The global dofs of the block of `V` at a point, and the process owning them.
1415 The dofs are in the global, unrolled numbering of the space `V` is an uncollapsed subspace of
1416 (or of `V`), as masters are given to :meth:`MultiPointConstraint.add_constraint` and
1417 :meth:`MultiPointConstraint.extend_masters`: one per component of `V`. Works for a space on a
1418 point mesh too.
1420 Args:
1421 V: The function space, or an uncollapsed subspace of it
1422 point: The point, of `gdim` or 3 coordinates
1423 distance_tol: The largest distance from the point to the dof. Defaults to `500` machine
1424 epsilon of the coordinate type of the mesh.
1426 Returns:
1427 The dofs, one per component, and the owning process. The same on every process.
1429 Raises:
1430 ValueError: On every process, unless exactly one block of dofs of `V` is at the point:
1431 none, or several, as in a discontinuous space or on a point mesh with coinciding
1432 points.
1434 Note:
1435 Collective.
1436 """
1437 mesh = V.mesh
1438 comm = mesh.comm
1439 gdim = mesh.geometry.dim
1440 if distance_tol is None:
1441 distance_tol = _tolerance(None, mesh.geometry.x.dtype)
1442 p = numpy.zeros(3, dtype=numpy.float64)
1443 given = numpy.asarray(point, dtype=numpy.float64).reshape(-1)
1444 p[: len(given)] = given
1446 def at_point(x):
1447 return numpy.linalg.norm(x[:gdim].T - p[:gdim], axis=1) <= distance_tol
1449 imap = V.dofmap.index_map
1450 parent_bs = V.dofmap.index_map_bs
1451 parent_dofs, components, _, bs = _marked_dofs(V, at_point)
1452 owned = parent_dofs < imap.size_local * parent_bs
1454 # The point must hold exactly one block of V, counted on the process owning it, which holds all
1455 # its components. A discontinuous space, or a point mesh with coinciding points, may have several.
1456 num_blocks = numpy.array([owned.sum() // bs], dtype=numpy.int64)
1457 comm.Allreduce(_MPI.IN_PLACE, num_blocks, op=_MPI.SUM)
1458 if num_blocks[0] != 1:
1459 raise ValueError(f"{num_blocks[0]} blocks of degrees of freedom of the space at {given}, not one")
1461 # The dofs and the owner, in one reduction
1462 found = numpy.full(bs + 1, -1, dtype=numpy.int64)
1463 found[components[owned]] = imap.local_range[0] * parent_bs + parent_dofs[owned]
1464 if owned.any():
1465 found[bs] = comm.rank
1466 comm.Allreduce(_MPI.IN_PLACE, found, op=_MPI.MAX)
1467 return found[:bs], int(found[bs])
1470def finalize_multipointconstraints(
1471 mpcs: Sequence[MultiPointConstraint], filter: Optional[numpy.floating] = None
1472) -> None:
1473 """
1474 Finalize the multi point constraints of several function spaces together.
1476 Entry ``k`` of ``mpcs`` constrains its own function space, for instance the ``k``-th block of a
1477 :class:`ufl.MixedFunctionSpace`. Each is finalized as by :meth:`MultiPointConstraint.finalize`,
1478 but the checks that need communication are reduced once for all of them, and the meshes of the
1479 spaces may be distinct, as long as they live on congruent communicators.
1481 Args:
1482 mpcs: The constraints to finalize. None may be finalized already, and they must all use the
1483 same ``dtype``.
1484 filter: See :meth:`MultiPointConstraint.finalize`. Applied to every constraint.
1486 Raises:
1487 ValueError: If the input is inconsistent, or if a dof is both a slave and constrained by a
1488 Dirichlet condition, a master is also a slave, or the meshes are on communicators
1489 of different size or rank order. Raised on every process.
1491 Note:
1492 Collective. Must be called by every process, with the constraints in the same order.
1493 """
1494 mpcs = list(mpcs)
1495 if len(mpcs) == 0:
1496 raise ValueError("At least one constraint is required")
1497 if len({id(mpc) for mpc in mpcs}) != len(mpcs):
1498 raise ValueError("The same constraint was given more than once")
1499 for mpc in mpcs:
1500 mpc._raise_if_finalized()
1501 dtype = numpy.dtype(mpcs[0]._dtype)
1502 if any(numpy.dtype(mpc._dtype) != dtype for mpc in mpcs):
1503 raise ValueError("All constraints must have the same dtype")
1504 if dtype.type not in (numpy.float32, numpy.float64, numpy.complex64, numpy.complex128):
1505 raise ValueError(f"Unsupported dtype {dtype} for coefficients")
1507 # An RBE3 constraint is built now, when the block of each space of its feet is known
1508 for mpc in mpcs:
1509 if len(mpc._rbe3) > 0:
1510 mpc._build_rbe3(mpcs)
1512 rhs_coeffs = []
1513 for mpc in mpcs:
1514 if mpc._rhs_coeffs is None:
1515 rhs_coeffs.append(numpy.zeros(0, dtype=dtype))
1516 else:
1517 num_dofs_local = mpc.V.dofmap.index_map_bs * (
1518 mpc.V.dofmap.index_map.size_local + mpc.V.dofmap.index_map.num_ghosts
1519 )
1520 rhs_coeffs.append(mpc._rhs_coeffs.x.array[:num_dofs_local].astype(dtype))
1522 def block_of(space: _fem.FunctionSpace) -> int:
1523 matches = [j for j, other in enumerate(mpcs) if other.V is space]
1524 if len(matches) == 0:
1525 matches = [j for j, other in enumerate(mpcs) if other.V.contains(space)]
1526 if len(matches) != 1:
1527 raise ValueError(
1528 "The master space of a constraint must be the function space, or an uncollapsed subspace of "
1529 "the function space, of exactly one of the constraints finalized together with it"
1530 )
1531 return matches[0]
1533 # The block of each master, from its code (see MultiPointConstraint.__init__). Every process
1534 # records the same spaces in the same order, so a space that is not one of the blocks raises
1535 # everywhere.
1536 master_blocks = []
1537 for k, mpc in enumerate(mpcs):
1538 codes = numpy.asarray(mpc._data.master_blocks)
1539 blocks = codes.copy()
1540 blocks[codes == -1] = k
1541 for s, space in enumerate(mpc._master_spaces):
1542 blocks[codes == -2 - s] = block_of(space)
1543 if len(mpc._master_spaces) == 0 and (codes == -1).all():
1544 blocks = numpy.zeros(0, dtype=numpy.int32)
1545 master_blocks.append(blocks)
1547 # Raises ValueError (as the C++ throws std::invalid_argument), identically on every process
1548 cpp_objects = dolfinx_mpc.cpp.mpc.create_multipointconstraints(
1549 [mpc.V._cpp_object for mpc in mpcs],
1550 # Copies: the arrays of the rows are read-only views of their C++ object
1551 [numpy.array(mpc._data.slaves) for mpc in mpcs],
1552 [numpy.array(mpc._data.masters) for mpc in mpcs],
1553 [mpc._data.coeffs.astype(dtype) for mpc in mpcs],
1554 [numpy.array(mpc._data.owners) for mpc in mpcs],
1555 [numpy.array(mpc._data.offsets) for mpc in mpcs],
1556 rhs_coeffs,
1557 [[bc._cpp_object for bc in mpc._bcs] for mpc in mpcs],
1558 master_blocks,
1559 filter,
1560 )
1562 # The block of each space on a spider mesh, matched before the spaces are replaced
1563 for mpc in mpcs:
1564 for entry in mpc._rbe2:
1565 entry[2] = next(j for j, other in enumerate(mpcs) if other.V is entry[0])
1567 for mpc, cpp_object in zip(mpcs, cpp_objects):
1568 mpc._cpp_object = cpp_object
1569 # Replace function space
1570 mpc.V = _fem.FunctionSpace(mpc.V.mesh, mpc.V.ufl_element(), cpp_object.function_space)
1571 mpc.finalized = True
1572 # Delete variables that are no longer required
1573 del (mpc._data, mpc._master_spaces)