Enhance IntegralForm: Non-symmetric mixed forms and None

CHANGELOG.md

@@ -10,6 +10,7 @@ All notable changes to this project will be documented in this file. The format
 - 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.
+- Add support for non-symmetric bilinear mixed forms in `IntegralForm`.
 
 ### Changed
 - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`.
@@ -19,6 +20,7 @@ All notable changes to this project will be documented in this file. The format
 - 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.
+- Integrate and assemble `None`-items in the `fun`-argument of `IntegralForm`. On integration, `None` is returned. `None` will be assembled to an emtpy sparse matrix.
 
 ### Fixed
 - Fix `Boundary(..., mode="and")` by ignoring any undefined axis.

src/felupe/assembly/_cartesian.py

@@ -104,7 +104,10 @@ class IntegralFormCartesian:
     """
 
     def __init__(self, fun, v, dV, u=None, grad_v=False, grad_u=False):
-        self.fun = np.ascontiguousarray(fun)
+        self.fun = fun
+        if self.fun is not None:
+            self.fun = np.ascontiguousarray(self.fun)
+
         self.dV = dV
 
         self.v = v
@@ -137,20 +140,22 @@ def assemble(self, values=None, parallel=False, out=None):
         if values is None:
             values = self.integrate(parallel=parallel, out=out)
 
-        permute = np.append(len(values.shape) - 1, range(len(values.shape) - 1)).astype(
-            int
-        )
+        if values is not None:
+            permute = np.append(
+                len(values.shape) - 1, range(len(values.shape) - 1)
+            ).astype(int)
 
-        # broadcast values of a uniform grid mesh
-        if values.size < self.indices[0].size:
-            new_shape = (*values.shape[:-1], self.v.region.mesh.ncells)
-            values = np.broadcast_to(values, new_shape)
+            # broadcast values of a uniform grid mesh
+            if values.size < self.indices[0].size:
+                new_shape = (*values.shape[:-1], self.v.region.mesh.ncells)
+                values = np.broadcast_to(values, new_shape)
 
-        res = sparsematrix(
-            (values.transpose(permute).ravel(), self.indices), shape=self.shape
-        )
+            return sparsematrix(
+                (values.transpose(permute).ravel(), self.indices), shape=self.shape
+            )
 
-        return res
+        else:
+            return sparsematrix(self.shape)
 
     def integrate(self, parallel=False, out=None):
         "Return evaluated (but not assembled) integrals."
@@ -160,6 +165,9 @@ def integrate(self, parallel=False, out=None):
         dV = self.dV
         fun = self.fun
 
+        if fun is None:
+            return None
+
         # plane strain
         # trim 3d vector-valued functions to the dimension of the field
         function_dimension = len(fun.shape) - 2

src/felupe/assembly/_integral.py

@@ -239,8 +239,8 @@ def __init__(self, fun, v, dV, u=None, grad_v=None, grad_u=None):
                 f = IntForm(fun=fun, v=v, dV=self.dV, grad_v=grad_v)
                 self.forms.append(f)
 
-        elif len(fun) == np.sum(1 + np.arange(self.nv)) and u is not None:
-            # BilinearForm
+        elif u is not None and len(fun) == np.sum(1 + np.arange(self.nv)):
+            # BilinearForm (symmetric, upper-triangle entries)
             self.mode = 2
             self.i, self.j = np.triu_indices(self.nv)
 
@@ -254,6 +254,23 @@ def __init__(self, fun, v, dV, u=None, grad_v=None, grad_u=None):
                     grad_u=self.grad_u[j],
                 )
                 self.forms.append(f)
+
+        elif u is not None and len(fun) == self.nv * self.nu:
+            # BilinearForm (non-symmetric)
+            self.mode = 3
+            self.i, self.j = [i.ravel() for i in np.indices((self.nv, self.nu))]
+
+            for a, (i, j) in enumerate(zip(self.i, self.j)):
+                f = IntForm(
+                    self.fun[a],
+                    v=self.v[i],
+                    dV=self.dV,
+                    u=self.u[j],
+                    grad_v=self.grad_v[i],
+                    grad_u=self.grad_u[j],
+                )
+                self.forms.append(f)
+
         else:
             raise ValueError("Unknown input format.")
 
@@ -267,14 +284,21 @@ def assemble(self, values=None, parallel=False, block=True, out=None):
             res.append(form.assemble(val, parallel=parallel, out=out))
 
         if block and self.mode == 2:
-            K = np.zeros((self.nv, self.nv), dtype=object)
+            K = np.zeros((self.nv, self.nu), dtype=object)
             for a, (i, j) in enumerate(zip(self.i, self.j)):
                 K[i, j] = res[a]
                 if i != j:
                     K[j, i] = res[a].T
 
             res = bmat(K).tocsr()
 
+        if block and self.mode == 3:
+            K = np.zeros((self.nv, self.nu), dtype=object)
+            for a, (i, j) in enumerate(zip(self.i, self.j)):
+                K[i, j] = res[a]
+
+            res = bmat(K).tocsr()
+
         if block and self.mode == 1:
             res = vstack(res).tocsr()
 

src/felupe/constitution/_mixed.py

@@ -214,8 +214,8 @@ def hessian(self, x, out=None):
         d2WdFdF = self.material.hessian([F, statevars], **kwargs)[0]
         d2WdFdF += p * detF * (dya(iFT, iFT) - cdya_il(iFT, iFT))
         d2WdFdp = detF * iFT
-        d2WdFdJ = np.zeros_like(F)
-        d2Wdpdp = np.zeros_like(p)
+        d2WdFdJ = None  # np.zeros_like(F)
+        d2Wdpdp = None  # np.zeros_like(p)
         d2WdpdJ = -np.ones_like(J)
         d2WdJdJ = self.d2UdJdJ(J, self.bulk) * np.ones_like(J)
         return [d2WdFdF, d2WdFdp, d2WdFdJ, d2Wdpdp, d2WdpdJ, d2WdJdJ]
@@ -623,7 +623,7 @@ def _hessian_pp(self, F, p, J):
         ndarray
             p,p - part of hessian
         """
-        return np.zeros_like(p)
+        return None  # np.zeros_like(p)
 
     def _hessian_JJ(self, F, p, J):
         """Linearization w.r.t. volume ratio of variation of