Add SolidBody(..., block=True, apply=None) (#920)

* Add `SolidBody(..., block=True, apply=None)`

to apply a final callable in `SolidBody.assemble.vector()` or `SolidBody.assemble.matrix()`.

* Add `SolidBody(..., multiplier=None)`

* Update test_mechanics.py

CHANGELOG.md

@@ -9,6 +9,7 @@ All notable changes to this project will be documented in this file. The format
 - Add a free-vibration modal analysis Step/Job `FreeVibration(items, boundaries)` with methods to evaluate `FreeVibration.evaluate()` and to extract `field, frequency = FreeVibration.extract(n)` its n-th result.
 - Add `Boundary.plot()` to plot the points and prescribed directions of a boundary.
 - Add `BoundaryDict` as a subclassed dict with methods to `plot()`, `screenshot()` and `imshow()`.
+- Add a new argument to apply a callable on the assembled vector/matrix of a solid body, `SolidBody(..., apply=None)`. This may be used to sum the list of sub-blocks instead of stacking them together, `SolidBody(..., block=False, apply=None)`. This is useful for mixed formulations where both the deformation gradient and the displacement values are required.
 
 ### Changed
 - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`.
@@ -17,6 +18,7 @@ All notable changes to this project will be documented in this file. The format
 - Enhance `Boundary` with added support for multiaxial prescribed values.
 - Enhance `math.linsteps(..., values=0)` with default values except for the column `axis` if `axis` is not None.
 - Link all field-values to the values of the first field if no other field is given in `FieldContainer.link()`.
+- Change the default arguments from `block=True` to `block=None` in `SolidBody.assemble.vector(block=None)` and `SolidBody.assemble.matrix(block=None)` and define `block` on creation, `SolidBody(..., block=True)` instead.
 
 ### Fixed
 - Fix `Boundary(..., mode="and")` by ignoring any undefined axis.

src/felupe/assembly/_integral.py

@@ -273,13 +273,12 @@ def assemble(self, values=None, parallel=False, block=True, out=None):
                 if i != j:
                     K[j, i] = res[a].T
 
-            return bmat(K).tocsr()
+            res = bmat(K).tocsr()
 
         if block and self.mode == 1:
-            return vstack(res).tocsr()
+            res = vstack(res).tocsr()
 
-        else:
-            return res
+        return res
 
     def integrate(self, parallel=False, out=None):
         if out is None:

src/felupe/mechanics/_solidbody.py

@@ -164,6 +164,13 @@ class SolidBody(Solid):
         Array of initial internal state variables (default is None).
     density : float or None, optional
         The density of the solid body.
+    block : bool, optional
+        Assemble a sparse matrix from sparse sub-blocks or assemble a sparse vector by
+        stacking sparse matrices vertically (row wise). Default is True.
+    apply : callable or None, optional
+        Apply a callable on the assembled vectors and sparse matrices. Default is None.
+    multiplier : float or None, optional
+        A scale factor for the assembled vector and matrix. Default is None.
 
     Notes
     -----
@@ -241,10 +248,21 @@ class SolidBody(Solid):
         methods for the assembly of sparse vectors/matrices.
     """
 
-    def __init__(self, umat, field, statevars=None, density=None):
+    def __init__(
+        self,
+        umat,
+        field,
+        statevars=None,
+        density=None,
+        block=True,
+        apply=None,
+        multiplier=None,
+    ):
         self.umat = umat
         self.field = field
         self.density = density
+        self.block = block
+        self.apply = apply
 
         self.results = Results(stress=True, elasticity=True)
         self.results.kinematics = self._extract(self.field)
@@ -264,7 +282,10 @@ def __init__(self, umat, field, statevars=None, density=None):
             )
 
         self.assemble = Assemble(
-            vector=self._vector, matrix=self._matrix, mass=self._mass
+            vector=self._vector,
+            matrix=self._matrix,
+            mass=self._mass,
+            multiplier=multiplier,
         )
 
         self.evaluate = Evaluate(
@@ -279,32 +300,61 @@ def __init__(self, umat, field, statevars=None, density=None):
         self._form = IntegralForm
 
     def _vector(
-        self, field=None, parallel=False, items=None, args=(), kwargs=None, block=True
+        self,
+        field=None,
+        parallel=False,
+        items=None,
+        args=(),
+        kwargs=None,
+        block=None,
+        apply=None,
     ):
         if kwargs is None:
             kwargs = {}
 
         if field is not None:
             self.field = field
 
+        if block is None:
+            block = self.block
+
+        if apply is None:
+            apply = self.apply
+
         self.results.stress = self._gradient(field, args=args, kwargs=kwargs)
         self.results.force = self._form(
             fun=self.results.stress[slice(items)],
             v=self.field,
             dV=self.field.region.dV,
         ).assemble(parallel=parallel, block=block)
 
+        if apply is not None:
+            self.results.force = apply(self.results.force)
+
         return self.results.force
 
     def _matrix(
-        self, field=None, parallel=False, items=None, args=(), kwargs=None, block=True
+        self,
+        field=None,
+        parallel=False,
+        items=None,
+        args=(),
+        kwargs=None,
+        block=None,
+        apply=None,
     ):
         if kwargs is None:
             kwargs = {}
 
         if field is not None:
             self.field = field
 
+        if block is None:
+            block = self.block
+
+        if apply is None:
+            apply = self.apply
+
         self.results.elasticity = self._hessian(field, args=args, kwargs=kwargs)
 
         form = self._form(
@@ -322,6 +372,9 @@ def _matrix(
             values=self.results.stiffness_values, block=block
         )
 
+        if apply is not None:
+            self.results.stiffness = apply(self.results.stiffness)
+
         return self.results.stiffness
 
     def _extract(self, field):

tests/test_mechanics.py

@@ -293,7 +293,7 @@ def test_solidbody_incompressible():
         assert np.allclose(t1, t2)
 
 
-def test_solidbody_axi():
+def test_solidbody_axi_incompressible():
     umat, u = pre_axi(bulk=None)
     b = fem.SolidBodyNearlyIncompressible(umat=umat, field=u, bulk=5000)
     b = fem.SolidBodyNearlyIncompressible(
@@ -338,7 +338,7 @@ def test_solidbody_axi():
         assert np.allclose(t1, t2)
 
 
-def test_solidbody_axi_incompressible():
+def test_solidbody_axi():
     umat, u = pre_axi()
     b = fem.SolidBody(umat=umat, field=u)
 
@@ -351,12 +351,18 @@ def test_solidbody_axi_incompressible():
         r2 = b.assemble.vector(**kwargs)
         assert np.allclose(r1.toarray(), r2.toarray())
 
+        r3 = b.assemble.vector(**kwargs, block=False, apply=sum)
+        assert np.allclose(r1.toarray(), r3.toarray())
+
         K1 = b.assemble.matrix(u, **kwargs)
         assert K1.shape == (18, 18)
 
         K2 = b.assemble.matrix(**kwargs)
         assert np.allclose(K1.toarray(), K2.toarray())
 
+        K3 = b.assemble.matrix(**kwargs, block=False, apply=sum)
+        assert np.allclose(K1.toarray(), K3.toarray())
+
         P1 = b.results.stress
         P2 = b.evaluate.gradient()
         P2 = b.evaluate.gradient(u)