Coverage for python/src/dolfinx_mpc/multipointconstraint.py: 82%
175 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-03 13:59 +0000
1# Copyright (C) 2020-2023 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, Tuple, Union
10from petsc4py import PETSc as _PETSc
12import dolfinx.cpp as _cpp
13import dolfinx.fem as _fem
14import dolfinx.mesh as _mesh
15import numpy
16import numpy.typing as npt
17from dolfinx import default_real_type, default_scalar_type
19import dolfinx_mpc.cpp
21from .dictcondition import create_dictionary_constraint
23_mpc_classes = Union[
24 dolfinx_mpc.cpp.mpc.MultiPointConstraint_double,
25 dolfinx_mpc.cpp.mpc.MultiPointConstraint_float,
26 dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_double,
27 dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_float,
28]
29_float_classes = Union[numpy.float32, numpy.float64, numpy.complex128, numpy.complex64]
30_float_array_types = Union[
31 npt.NDArray[numpy.float32],
32 npt.NDArray[numpy.float64],
33 npt.NDArray[numpy.complex64],
34 npt.NDArray[numpy.complex128],
35]
36_mpc_data_classes = Union[
37 dolfinx_mpc.cpp.mpc.mpc_data_double,
38 dolfinx_mpc.cpp.mpc.mpc_data_float,
39 dolfinx_mpc.cpp.mpc.mpc_data_complex_double,
40 dolfinx_mpc.cpp.mpc.mpc_data_complex_float,
41]
44class MPCData:
45 _cpp_object: _mpc_data_classes
47 def __init__(
48 self,
49 slaves: npt.NDArray[numpy.int32],
50 masters: npt.NDArray[numpy.int64],
51 coeffs: _float_array_types,
52 owners: npt.NDArray[numpy.int32],
53 offsets: npt.NDArray[numpy.int32],
54 ):
55 if coeffs.dtype.type == numpy.float32:
56 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_float(slaves, masters, coeffs, owners, offsets)
57 elif coeffs.dtype.type == numpy.float64:
58 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_double(slaves, masters, coeffs, owners, offsets)
59 elif coeffs.dtype.type == numpy.complex64:
60 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_complex_float(slaves, masters, coeffs, owners, offsets)
61 elif coeffs.dtype.type == numpy.complex128:
62 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_complex_double(slaves, masters, coeffs, owners, offsets)
63 else:
64 raise ValueError("Unsupported dtype {coeffs.dtype.type} for coefficients")
66 @property
67 def slaves(self):
68 return self._cpp_object.slaves
70 @property
71 def masters(self):
72 return self._cpp_object.masters
74 @property
75 def coeffs(self):
76 return self._cpp_object.coeffs
78 @property
79 def owners(self):
80 return self._cpp_object.owners
82 @property
83 def offsets(self):
84 return self._cpp_object.offsets
87class MultiPointConstraint:
88 """
89 Hold data for multi point constraint relation ships,
90 including new index maps for local assembly of matrices and vectors.
92 Args:
93 V: The function space
94 dtype: The dtype of the underlying functions
95 """
97 _slaves: npt.NDArray[numpy.int32]
98 _masters: npt.NDArray[numpy.int64]
99 _coeffs: _float_array_types
100 _owners: npt.NDArray[numpy.int32]
101 _offsets: npt.NDArray[numpy.int32]
102 V: _fem.FunctionSpace
103 finalized: bool
104 _cpp_object: _mpc_classes
105 _dtype: numpy.floating
106 __slots__ = tuple(__annotations__)
108 def __init__(self, V: _fem.FunctionSpace, dtype: numpy.floating = default_scalar_type):
109 self._slaves = numpy.array([], dtype=numpy.int32)
110 self._masters = numpy.array([], dtype=numpy.int64)
111 self._coeffs = numpy.array([], dtype=dtype) # type: ignore
112 self._owners = numpy.array([], dtype=numpy.int32)
113 self._offsets = numpy.array([0], dtype=numpy.int32)
114 self.V = V
115 self.finalized = False
116 self._dtype = dtype
118 def add_constraint(
119 self,
120 V: _fem.FunctionSpace,
121 slaves: npt.NDArray[numpy.int32],
122 masters: npt.NDArray[numpy.int64],
123 coeffs: _float_array_types,
124 owners: npt.NDArray[numpy.int32],
125 offsets: npt.NDArray[numpy.int32],
126 ):
127 """
128 Add new constraint given by numpy arrays.
130 Args:
131 V: The function space for the constraint
132 slaves: List of all slave dofs (using local dof numbering) on this process
133 masters: List of all master dofs (using global dof numbering) on this process
134 coeffs: The coefficients corresponding to each master.
135 owners: The process each master is owned by.
136 offsets: Array indicating the location in the masters array for the i-th slave
137 in the slaves arrays, i.e.
139 .. highlight:: python
140 .. code-block:: python
142 masters_of_owned_slave[i] = masters[offsets[i]:offsets[i+1]]
144 """
145 assert V == self.V
146 self._already_finalized()
148 if len(slaves) > 0:
149 self._offsets = numpy.append(self._offsets, offsets[1:] + len(self._masters))
150 self._slaves = numpy.append(self._slaves, slaves)
151 self._masters = numpy.append(self._masters, masters)
152 self._coeffs = numpy.append(self._coeffs, coeffs)
153 self._owners = numpy.append(self._owners, owners)
155 def add_constraint_from_mpc_data(self, V: _fem.FunctionSpace, mpc_data: Union[_mpc_data_classes, MPCData]):
156 """
157 Add new constraint given by an `dolfinc_mpc.cpp.mpc.mpc_data`-object
158 """
159 self._already_finalized()
160 self.add_constraint(
161 V,
162 mpc_data.slaves,
163 mpc_data.masters,
164 mpc_data.coeffs,
165 mpc_data.owners,
166 mpc_data.offsets,
167 )
169 def finalize(self) -> None:
170 """
171 Finializes the multi point constraint. After this function is called, no new constraints can be added
172 to the constraint. This function creates a map from the cells (local to index) to the slave degrees of
173 freedom and builds a new index map and function space where unghosted master dofs are added as ghosts.
174 """
175 self._already_finalized()
176 self._coeffs.astype(numpy.dtype(self._dtype))
177 # Initialize C++ object and create slave->cell maps
178 if self._dtype == numpy.float32:
179 self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_float(
180 self.V._cpp_object,
181 self._slaves,
182 self._masters,
183 self._coeffs.astype(self._dtype),
184 self._owners,
185 self._offsets,
186 )
187 elif self._dtype == numpy.float64:
188 self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_double(
189 self.V._cpp_object,
190 self._slaves,
191 self._masters,
192 self._coeffs.astype(self._dtype),
193 self._owners,
194 self._offsets,
195 )
196 elif self._dtype == numpy.complex64:
197 self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_float(
198 self.V._cpp_object,
199 self._slaves,
200 self._masters,
201 self._coeffs.astype(self._dtype),
202 self._owners,
203 self._offsets,
204 )
206 elif self._dtype == numpy.complex128:
207 self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_double(
208 self.V._cpp_object,
209 self._slaves,
210 self._masters,
211 self._coeffs.astype(self._dtype),
212 self._owners,
213 self._offsets,
214 )
215 else:
216 raise ValueError("Unsupported dtype {coeffs.dtype.type} for coefficients")
218 # Replace function space
219 self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
221 self.finalized = True
222 # Delete variables that are no longer required
223 del (self._slaves, self._masters, self._coeffs, self._owners, self._offsets)
225 def create_periodic_constraint_topological(
226 self,
227 V: _fem.FunctionSpace,
228 meshtag: _mesh.MeshTags,
229 tag: int,
230 relation: Callable[[numpy.ndarray], numpy.ndarray],
231 bcs: List[_fem.DirichletBC],
232 scale: _float_classes = default_scalar_type(1.0),
233 tol: _float_classes = 500 * numpy.finfo(default_real_type).eps,
234 ):
235 """
236 Create periodic condition for all closure dofs of on all entities in `meshtag` with value `tag`.
237 :math:`u(x_i) = scale * u(relation(x_i))` for all of :math:`x_i` on marked entities.
239 Args:
240 V: The function space to assign the condition to. Should either be the space of the MPC or a sub space.
241 meshtag: MeshTag for entity to apply the periodic condition on
242 tag: Tag indicating which entities should be slaves
243 relation: Lambda-function describing the geometrical relation
244 bcs: Dirichlet boundary conditions for the problem (Periodic constraints will be ignored for these dofs)
245 scale: Float for scaling bc
246 tol: Tolerance for adding scaled basis values to MPC. Any contribution that is less than this value
247 is ignored. The tolerance is also added as padding for the bounding box trees and corresponding
248 collision searches to determine periodic degrees of freedom.
249 """
250 bcs_ = [bc._cpp_object for bc in bcs]
251 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
252 scale = scale.item() # type: ignore
253 if V is self.V:
254 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological(
255 self.V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, False, float(tol)
256 )
257 elif self.V.contains(V):
258 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological(
259 V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, True, float(tol)
260 )
261 else:
262 raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")
263 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
265 def create_periodic_constraint_geometrical(
266 self,
267 V: _fem.FunctionSpace,
268 indicator: Callable[[numpy.ndarray], numpy.ndarray],
269 relation: Callable[[numpy.ndarray], numpy.ndarray],
270 bcs: List[_fem.DirichletBC],
271 scale: _float_classes = default_scalar_type(1.0),
272 tol: _float_classes = 500 * numpy.finfo(default_real_type).eps,
273 ):
274 """
275 Create a periodic condition for all degrees of freedom whose physical location satisfies
276 :math:`indicator(x_i)==True`, i.e.
277 :math:`u(x_i) = scale * u(relation(x_i))` for all :math:`x_i`
279 Args:
280 V: The function space to assign the condition to. Should either be the space of the MPC or a sub space.
281 indicator: Lambda-function to locate degrees of freedom that should be slaves
282 relation: Lambda-function describing the geometrical relation to master dofs
283 bcs: Dirichlet boundary conditions for the problem
284 (Periodic constraints will be ignored for these dofs)
285 scale: Float for scaling bc
286 tol: Tolerance for adding scaled basis values to MPC. Any contribution that is less than this value
287 is ignored. The tolerance is also added as padding for the bounding box trees and corresponding
288 collision searches to determine periodic degrees of freedom.
289 """
290 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
291 scale = scale.item() # type: ignore
292 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
293 if V is self.V:
294 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical(
295 self.V._cpp_object, indicator, relation, bcs, scale, False, float(tol)
296 )
297 elif self.V.contains(V):
298 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical(
299 V._cpp_object, indicator, relation, bcs, scale, True, float(tol)
300 )
301 else:
302 raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")
303 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
305 def create_slip_constraint(
306 self,
307 space: _fem.FunctionSpace,
308 facet_marker: Tuple[_mesh.MeshTags, int],
309 v: _fem.Function,
310 bcs: List[_fem.DirichletBC] = [],
311 ):
312 """
313 Create a slip constraint :math:`u \\cdot v=0` over the entities defined in `facet_marker` with the given index.
315 Args:
316 space: Function space (possible sub space) for the current constraint
317 facet_marker: Tuple containomg the mesh tag and marker used to locate degrees of freedom
318 v: Function containing the directional vector to dot your slip condition (most commonly a normal vector)
319 bcs: List of Dirichlet BCs (slip conditions will be ignored on these dofs)
321 Examples:
322 Create constaint :math:`u\\cdot n=0` of all indices in `mt` marked with `i`
324 .. highlight:: python
325 .. code-block:: python
327 V = dolfinx.fem.functionspace(mesh, ("CG", 1))
328 mpc = MultiPointConstraint(V)
329 n = dolfinx.fem.Function(V)
330 mpc.create_slip_constaint(V, (mt, i), n)
332 Create slip constaint for a mixed function space:
334 .. highlight:: python
335 .. code-block:: python
337 cellname = mesh.ufl_cell().cellname()
338 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
339 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
340 me = basix.ufl.mixed_element([Ve, Qe])
341 W = dolfinx.fem.functionspace(mesh, me)
342 mpc = MultiPointConstraint(W)
343 n_space, _ = W.sub(0).collapse()
344 normal = dolfinx.fem.Function(n_space)
345 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[])
347 A slip condition cannot be applied on the same degrees of freedom as a Dirichlet BC, and therefore
348 any Dirichlet bc for the space of the multi point constraint should be supplied.
350 .. highlight:: python
351 .. code-block:: python
353 cellname = mesh.ufl_cell().cellname()
354 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
355 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
356 me = basix.ufl.mixed_element([Ve, Qe])
357 W = dolfinx.fem.functionspace(mesh, me)
358 mpc = MultiPointConstraint(W)
359 n_space, _ = W.sub(0).collapse()
360 normal = Function(n_space)
361 bc = dolfinx.fem.dirichletbc(inlet_velocity, dofs, W.sub(0))
362 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[bc])
363 """
364 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
365 if space is self.V:
366 sub_space = False
367 elif self.V.contains(space):
368 sub_space = True
369 else:
370 raise ValueError("Input space has to be a sub space of the MPC space")
371 mpc_data = dolfinx_mpc.cpp.mpc.create_slip_condition(
372 space._cpp_object,
373 facet_marker[0]._cpp_object,
374 facet_marker[1],
375 v._cpp_object,
376 bcs,
377 sub_space,
378 )
379 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
381 def create_general_constraint(
382 self,
383 slave_master_dict: Dict[bytes, Dict[bytes, float]],
384 subspace_slave: Optional[int] = None,
385 subspace_master: Optional[int] = None,
386 ):
387 """
388 Args:
389 V: The function space
390 slave_master_dict: Nested dictionary, where the first key is the bit representing the slave dof's
391 coordinate in the mesh. The item of this key is a dictionary, where each key of this dictionary
392 is the bit representation of the master dof's coordinate, and the item the coefficient for
393 the MPC equation.
394 subspace_slave: If using mixed or vector space, and only want to use dofs from a sub space
395 as slave add index here
396 subspace_master: Subspace index for mixed or vector spaces
398 Example:
399 If the dof `D` located at `[d0, d1]` should be constrained to the dofs
400 `E` and `F` at `[e0, e1]` and `[f0, f1]` as :math:`D = \\alpha E + \\beta F`
401 the dictionary should be:
403 .. highlight:: python
404 .. code-block:: python
406 {numpy.array([d0, d1], dtype=mesh.geometry.x.dtype).tobytes():
407 {numpy.array([e0, e1], dtype=mesh.geometry.x.dtype).tobytes(): alpha,
408 numpy.array([f0, f1], dtype=mesh.geometry.x.dtype).tobytes(): beta}}
409 """
410 slaves, masters, coeffs, owners, offsets = create_dictionary_constraint(
411 self.V, slave_master_dict, subspace_slave, subspace_master
412 )
413 self.add_constraint(self.V, slaves, masters, coeffs, owners, offsets)
415 def create_contact_slip_condition(
416 self,
417 meshtags: _mesh.MeshTags,
418 slave_marker: int,
419 master_marker: int,
420 normal: _fem.Function,
421 eps2: float = 1e-20,
422 ):
423 """
424 Create a slip condition between two sets of facets marker with individual markers.
425 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
426 The condition created is :math:`u_s \\cdot normal_s = u_m \\cdot normal_m` where `s` is the
427 restriction to the slave facets, `m` to the master facets.
429 Args:
430 meshtags: The meshtags of the set of facets to tie together
431 slave_marker: The marker of the slave facets
432 master_marker: The marker of the master facets
433 normal: The function used in the dot-product of the constraint
434 eps2: The tolerance for the squared distance between cells to be considered as a collision
435 """
436 if isinstance(eps2, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
437 eps2 = eps2.item() # type: ignore
438 mpc_data = dolfinx_mpc.cpp.mpc.create_contact_slip_condition(
439 self.V._cpp_object,
440 meshtags._cpp_object,
441 slave_marker,
442 master_marker,
443 normal._cpp_object,
444 eps2,
445 )
446 self.add_constraint_from_mpc_data(self.V, mpc_data)
448 def create_contact_inelastic_condition(
449 self,
450 meshtags: _cpp.mesh.MeshTags_int32,
451 slave_marker: int,
452 master_marker: int,
453 eps2: float = 1e-20,
454 ):
455 """
456 Create a contact inelastic condition between two sets of facets marker with individual markers.
457 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
458 The condition created is :math:`u_s = u_m` where `s` is the restriction to the
459 slave facets, `m` to the master facets.
461 Args:
462 meshtags: The meshtags of the set of facets to tie together
463 slave_marker: The marker of the slave facets
464 master_marker: The marker of the master facets
465 eps2: The tolerance for the squared distance between cells to be considered as a collision
466 """
467 if isinstance(eps2, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
468 eps2 = eps2.item() # type: ignore
469 mpc_data = dolfinx_mpc.cpp.mpc.create_contact_inelastic_condition(
470 self.V._cpp_object, meshtags._cpp_object, slave_marker, master_marker, eps2
471 )
472 self.add_constraint_from_mpc_data(self.V, mpc_data)
474 @property
475 def is_slave(self) -> numpy.ndarray:
476 """
477 Returns a vector of integers where the ith entry indicates if a degree of freedom (local to process) is a slave.
478 """
479 self._not_finalized()
480 return self._cpp_object.is_slave
482 @property
483 def slaves(self):
484 """
485 Returns the degrees of freedom for all slaves local to process
486 """
487 self._not_finalized()
488 return self._cpp_object.slaves
490 @property
491 def masters(self) -> _cpp.graph.AdjacencyList_int32:
492 """
493 Returns an adjacency-list whose ith node corresponds to
494 a degree of freedom (local to process), and links the corresponding master dofs (local to process).
496 Examples:
498 .. highlight:: python
499 .. code-block:: python
501 masters = mpc.masters
502 masters_of_dof_i = masters.links(i)
503 """
504 self._not_finalized()
505 return self._cpp_object.masters
507 def coefficients(self) -> _float_array_types:
508 """
509 Returns a vector containing the coefficients for the constraint, and the corresponding offsets
510 for the ith degree of freedom.
512 Examples:
514 .. highlight:: python
515 .. code-block:: python
517 coeffs, offsets = mpc.coefficients()
518 coeffs_of_slave_i = coeffs[offsets[i]:offsets[i+1]]
519 """
520 self._not_finalized()
521 return self._cpp_object.coefficients()
523 @property
524 def num_local_slaves(self):
525 """
526 Return the number of slaves owned by the current process.
527 """
528 self._not_finalized()
529 return self._cpp_object.num_local_slaves
531 @property
532 def cell_to_slaves(self):
533 """
534 Returns an `dolfinx.cpp.graph.AdjacencyList_int32` whose ith node corresponds to
535 the ith cell (local to process), and links the corresponding slave degrees of
536 freedom in the cell (local to process).
538 Examples:
540 .. highlight:: python
541 .. code-block:: python
543 cell_to_slaves = mpc.cell_to_slaves()
544 slaves_in_cell_i = cell_to_slaves.links(i)
545 """
546 self._not_finalized()
547 return self._cpp_object.cell_to_slaves
549 @property
550 def function_space(self):
551 """
552 Return the function space for the multi-point constraint with the updated index map
553 """
554 self._not_finalized()
555 return self.V
557 def backsubstitution(self, u: Union[_fem.Function, _PETSc.Vec]) -> None: # type: ignore
558 """
559 For a Function, impose the multi-point constraint by backsubstiution.
560 This function is used after solving the reduced problem to obtain the values
561 at the slave degrees of freedom
563 .. note::
564 It is the users responsibility to destroy the PETSc vector
566 Args:
567 u: The input function
568 """
569 try:
570 self._cpp_object.backsubstitution(u.x.array)
571 u.x.scatter_forward()
572 except AttributeError:
573 with u.localForm() as vector_local:
574 self._cpp_object.backsubstitution(vector_local.array_w)
575 u.ghostUpdate(addv=_PETSc.InsertMode.INSERT, mode=_PETSc.ScatterMode.FORWARD) # type: ignore
577 def homogenize(self, u: _fem.Function) -> None:
578 """
579 For a vector, homogenize (set to zero) the vector components at the multi-point
580 constraint slave DoF indices. This is particularly useful for nonlinear problems.
582 Args:
583 u: The input vector
584 """
585 self._cpp_object.homogenize(u.x.array)
586 u.x.scatter_forward()
588 def _already_finalized(self):
589 """
590 Check if we have already finalized the multi point constraint
591 """
592 if self.finalized:
593 raise RuntimeError("MultiPointConstraint has already been finalized")
595 def _not_finalized(self):
596 """
597 Check if we have finalized the multi point constraint
598 """
599 if not self.finalized:
600 raise RuntimeError("MultiPointConstraint has not been finalized")