adtzlr · adtzlr · Aug 7, 2024 · Aug 7, 2024 · Aug 7, 2024 · Aug 7, 2024
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -12,9 +12,12 @@ All notable changes to this project will be documented in this file. The format
 - Add optional keyword-argument `SolidBody.assemble.matrix(block=True)`, also for `SolidBody.assemble.vector(block=True)` and the assemble-methods of `SolidBodyNearlyIncompressible`. If `block=False`, these methods will assemble a list of the upper-triangle sub block-vectors/-matrices instead of the combined block-vector/-matrix.
 - Add `SolidBodyForce` as a replacement for `SolidBodyGravity` with more general keyword arguments (`"values"` instead of `"gravity"`, `"scale"` instead of `"density"`).
 - Add `Laplace` to be used as user-material in a solid body. Works with scalar- and vector-valued fields.
+- Add an optional `Assemble(..., multiplier=None)` argument which is used in external items like  `SolidBodyForce`, `SolidBodyGravity` and `PointLoad` and is applied in `newtonrhapson(items, ...)`.
 
 ### Changed
 - Change the internal initialization of `field = Field(region, values=1, dtype=None)` values from `field.values = np.ones(shape) * values` to `field = np.full(shape, fill_value=values, dtype=dtype)`. This enforces `field = Field(region, values=1)` to return the gradient array with data-type `int` which was of type `float` before.
+- Initialize empty matrices of `SolidBodyForce`, `SolidBodyGravity` and `PointLoad` with `dtype=float`.
+- Don't multiply the assembled vectors of `SolidBodyForce`, `SolidBodyGravity` and `PointLoad` by `-1.0`.
 
 ### Deprecated
 - Deprecate `SolidBodyGravity`, `SolidBodyForce` should be used instead.

diff --git a/docs/howto/solid.rst b/docs/howto/solid.rst
@@ -8,22 +8,23 @@ Solid Body
 
 The generation of internal force vectors and stiffness matrices of solid bodies are provided as assembly-methods of a :class:`~felupe.SolidBody` or a :class:`~felupe.SolidBodyNearlyIncompressible`.
 
-..  code-block:: python
+..  pyvista-plot::
+    :context:
 
     import felupe as fem
 
     mesh = fem.Cube(n=6)
     region = fem.RegionHexahedron(mesh)
     field = fem.FieldContainer([fem.Field(region, dim=3)])
-    
+
     # a solid body
     body = fem.SolidBody(umat=fem.NeoHooke(mu=1, bulk=2), field=field)
-    
+
     # a (nearly) incompressible solid body (to be used with quads and hexahedrons)
     body = fem.SolidBodyNearlyIncompressible(umat=fem.NeoHooke(mu=1), field=field, bulk=5000)
-    
-    internal_force = body.assemble.vector(field, parallel=False, jit=False)
-    stiffness_matrix = body.assemble.matrix(field, parallel=False, jit=False)
+
+    internal_force = body.assemble.vector(field, parallel=False)
+    stiffness_matrix = body.assemble.matrix(field, parallel=False)
 
 
 During assembly, several results are stored, e.g. the gradient of the strain energy density function per unit undeformed volume w.r.t. the deformation gradient (first Piola-Kirchhoff stress tensor). Other results are the deformation gradient or the fourth-order elasticity tensor associated to the first Piola-Kirchhoff stress tensor.
@@ -36,7 +37,8 @@ During assembly, several results are stored, e.g. the gradient of the strain ene
 
     \mathbb{A} &= \frac{\partial^2 \psi(\boldsymbol{C}(\boldsymbol{F}))}{\partial \boldsymbol{F}\ \partial \boldsymbol{F}}
 
-..  code-block:: python
+..  pyvista-plot::
+    :context:
 
     F = body.results.kinematics
     P = body.results.stress
@@ -54,8 +56,9 @@ The Cauchy stress tensor, as well as the gradient and the hessian of the strain
     \boldsymbol{\sigma} &= \frac{1}{J} \boldsymbol{P} \boldsymbol{F}^T
 
 
-..  code-block:: python
-
+..  pyvista-plot::
+    :context:
+
     P = body.evaluate.gradient(field)
     A = body.evaluate.hessian(field)
     s = body.evaluate.cauchy_stress(field)
@@ -72,11 +75,12 @@ The generation of external force vectors or stiffness matrices of body forces ac
     \delta W_{ext} = \int_V \delta \boldsymbol{u} \cdot \rho \boldsymbol{g} \ dV
 
 
-..  code-block:: python
-
+..  pyvista-plot::
+    :context:
+
     body = fem.SolidBodyForce(field=field, values=[9810, 0, 0], scale=7.85e-9)
-    
-    force_gravity = body.assemble.vector(field, parallel=False, jit=False)
+
+    force_gravity = body.assemble.vector(field, parallel=False)
 
 
 Pressure Boundary on a Solid Body
@@ -86,44 +90,44 @@ The generation of force vectors or stiffness matrices of pressure boundaries on
 
 ..  math::
 
-    \delta W_{ext} = \int_{\partial V} \delta \boldsymbol{u} \cdot p \ J \boldsymbol{F}^{-T} \ d\boldsymbol{A}
+    \delta W_{ext} = \int_{\partial V} \delta \boldsymbol{u} \cdot (-p) \ J \boldsymbol{F}^{-T} \ d\boldsymbol{A}
+
+
+..  pyvista-plot::
+    :context:
 
-..  code-block:: python
-
     region_pressure = fem.RegionHexahedronBoundary(
         mesh=mesh,
         only_surface=True, # select only faces on the outline
         mask=mesh.points[:, 0] == 0, # select a subset of faces on the surface
     )
-    
+
     field_boundary = fem.FieldContainer([fem.Field(region_pressure, dim=3)])
     field_boundary.link(field)
-    
+
     body_pressure = fem.SolidBodyPressure(field=field_boundary)
-
-    force_pressure = body_pressure.assemble.vector(
-        field=field_boundary, parallel=False, jit=False
-    )
-
+
+    force_pressure = body_pressure.assemble.vector(field=field_boundary, parallel=False)
     stiffness_matrix_pressure = body_pressure.assemble.matrix(
-        field=field_boundary, parallel=False, jit=False
+        field=field_boundary, parallel=False
     )
 
 
 For axisymmetric problems the boundary region has to be created with the attribute ``ensure_3d=True``.
 
-..  code-block:: python
-
+..  pyvista-plot::
+    :context:
+
     mesh = fem.Rectangle(a=(0, 30), b=(20, 40), n=(21, 11))
     region = fem.RegionQuad(mesh)
-    
+
     region_pressure = fem.RegionQuadBoundary(
         mesh=mesh,
         only_surface=True, # select only faces on the outline
         mask=mesh.points[:, 0] == 0, # select a subset of faces on the surface
         ensure_3d=True, # flag for axisymmetric boundary region
     )
-    
+
     field = fem.FieldContainer([fem.FieldAxisymmetric(region)])
     field_boundary = fem.FieldContainer([fem.FieldAxisymmetric(region_pressure)])
     field_boundary.link(field)
diff --git a/docs/howto/solvers.rst b/docs/howto/solvers.rst
@@ -21,7 +21,9 @@ Solvers from SciPy Sparse:
 
 .. tab:: SciPy Sparse (iterative)
 
-   **Note**: ``minres`` may be replaced by another iterative method.
+   ..  note::
+
+       ``minres`` may be replaced by another iterative method.
 
    ..  code-block:: python
 
@@ -45,10 +47,6 @@ Solvers from external packages:
    ..  code-block:: python
 
        from pypardiso import spsolve as solver
-
-       # undocumented, untested workaround if multiple blas libaries are installed
-       # import os
-       # os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
 
 .. tab:: PyPardiso (direct, symmetric)
 

diff --git a/src/felupe/field/_base.py b/src/felupe/field/_base.py
@@ -179,7 +179,7 @@ def grad(self, sym=False, dtype=None, out=None):
             return g
 
     def interpolate(self, dtype=None, out=None):
-        """Interpolate field values located at mesh-points to the quadrature points
+        r"""Interpolate field values located at mesh-points to the quadrature points
         ``q`` of cells ``c`` in the region.
 
         ..  math::

diff --git a/src/felupe/mechanics/_helpers.py b/src/felupe/mechanics/_helpers.py
@@ -24,11 +24,12 @@
 
 
 class Assemble:
-    "A class with assembly methods of a SolidBody."
+    "A class with methods for assembling vectors and matrices of a SolidBody."
 
-    def __init__(self, vector, matrix):
+    def __init__(self, vector, matrix, multiplier=None):
         self.vector = vector
         self.matrix = matrix
+        self.multiplier = multiplier
 
 
 class Evaluate:

diff --git a/src/felupe/mechanics/_pointload.py b/src/felupe/mechanics/_pointload.py
@@ -40,13 +40,6 @@ class PointLoad:
         A flag to multiply the assembled vector and matrix by a scaling factor of
         :math:`2 \pi` (default is False).
 
-    Notes
-    -----
-    .. warning::
-
-       The assembled vector is returned with a negative sign because this is considered
-       as an external quantity.
-
     Examples
     --------
     ..  pyvista-plot::
@@ -65,9 +58,9 @@ class PointLoad:
         >>>
         >>> vector = load.assemble.vector()
         >>> vector.toarray()
-        array([[ 0.],
-               [-3.],
-               [-5.]])
+        array([[0.],
+               [3.],
+               [5.]])
     """
 
     def __init__(self, field, points, values=None, apply_on=0, axisymmetric=False):
@@ -83,7 +76,9 @@ def __init__(self, field, points, values=None, apply_on=0, axisymmetric=False):
         self.axisymmetric = axisymmetric
 
         self.results = Results()
-        self.assemble = Assemble(vector=self._vector, matrix=self._matrix)
+        self.assemble = Assemble(
+            vector=self._vector, matrix=self._matrix, multiplier=-1.0
+        )
 
     def update(self, values):
         self.__init__(self.field, self.points, values, self.apply_on, self.axisymmetric)
@@ -104,13 +99,13 @@ def _vector(self, field=None, parallel=False):
             np.concatenate([f.ravel() for f in force]).reshape(-1, 1)
         )
 
-        return -self.results.force
+        return self.results.force
 
     def _matrix(self, field=None, parallel=False):
         if field is not None:
             self.field = field
 
         n = np.sum(self.field.fieldsizes)
-        self.results.stiffness = csr_matrix(([0], ([0], [0])), shape=(n, n))
+        self.results.stiffness = csr_matrix(([0.0], ([0], [0])), shape=(n, n))
 
         return self.results.stiffness
diff --git a/src/felupe/mechanics/_solidbody_force.py b/src/felupe/mechanics/_solidbody_force.py
@@ -76,7 +76,9 @@ class SolidBodyForce:
     def __init__(self, field, values=None, scale=1.0):
         self.field = field
         self.results = Results(stress=False, elasticity=False)
-        self.assemble = Assemble(vector=self._vector, matrix=self._matrix)
+        self.assemble = Assemble(
+            vector=self._vector, matrix=self._matrix, multiplier=-1.0
+        )
         self._form = IntegralForm
 
         self.results.values = np.zeros(self.field[0].dim)
@@ -106,13 +108,13 @@ def _vector(self, field=None, parallel=False):
         if len(self.field) > 1:
             self.results.force.resize(np.sum(self.field.fieldsizes), 1)
 
-        return -self.results.force
+        return self.results.force
 
     def _matrix(self, field=None, parallel=False):
         if field is not None:
             self.field = field
 
         n = np.sum(self.field.fieldsizes)
-        self.results.stiffness = csr_matrix(([0], ([0], [0])), shape=(n, n))
+        self.results.stiffness = csr_matrix(([0.0], ([0], [0])), shape=(n, n))
 
         return self.results.stiffness
diff --git a/src/felupe/mechanics/_solidbody_gravity.py b/src/felupe/mechanics/_solidbody_gravity.py
@@ -83,7 +83,9 @@ def __init__(self, field, gravity=None, density=1.0):
 
         self.field = field
         self.results = Results(stress=False, elasticity=False)
-        self.assemble = Assemble(vector=self._vector, matrix=self._matrix)
+        self.assemble = Assemble(
+            vector=self._vector, matrix=self._matrix, multiplier=-1.0
+        )
         self._form = IntegralForm
 
         self.results.gravity = np.zeros(self.field[0].dim)
@@ -113,13 +115,13 @@ def _vector(self, field=None, parallel=False):
         if len(self.field) > 1:
             self.results.force.resize(np.sum(self.field.fieldsizes), 1)
 
-        return -self.results.force
+        return self.results.force
 
     def _matrix(self, field=None, parallel=False):
         if field is not None:
             self.field = field
 
         n = np.sum(self.field.fieldsizes)
-        self.results.stiffness = csr_matrix(([0], ([0], [0])), shape=(n, n))
+        self.results.stiffness = csr_matrix(([0.0], ([0], [0])), shape=(n, n))
 
         return self.results.stiffness
diff --git a/src/felupe/mechanics/_solidbody_pressure.py b/src/felupe/mechanics/_solidbody_pressure.py
@@ -41,7 +41,7 @@ class SolidBodyPressure:
     ..  math::
 
         \delta W_{ext} = \int_{\partial V}
-            \delta \boldsymbol{u} \cdot p J \boldsymbol{F}^{-T} \ d\boldsymbol{A}
+            \delta \boldsymbol{u} \cdot (-p) J \boldsymbol{F}^{-T} \ d\boldsymbol{A}
 
     Examples
     --------
@@ -87,7 +87,9 @@ def __init__(self, field, pressure=None):
         if pressure is not None:
             self.results.pressure = pressure
 
-        self.assemble = Assemble(vector=self._vector, matrix=self._matrix)
+        self.assemble = Assemble(
+            vector=self._vector, matrix=self._matrix, multiplier=-1.0
+        )
         self._area_change = AreaChange()
 
     def update(self, pressure):
@@ -113,7 +115,7 @@ def _vector(self, field=None, pressure=None, parallel=False, resize=None):
         if pressure is not None:
             self.results.pressure = pressure
 
-        fun[0] *= self.results.pressure
+        fun[0] *= -self.results.pressure
 
         self.results.force = IntegralForm(
             fun=fun, v=self.field, dV=self.field.region.dV, grad_v=[False]
@@ -138,7 +140,7 @@ def _matrix(self, field=None, pressure=None, parallel=False, resize=None):
         if pressure is not None:
             self.results.pressure = pressure
 
-        fun[0] *= self.results.pressure
+        fun[0] *= -self.results.pressure
 
         self.results.stiffness = IntegralForm(
             fun=fun,

diff --git a/src/felupe/mesh/_mesh.py b/src/felupe/mesh/_mesh.py
@@ -1064,7 +1064,7 @@ def fill_between(self, other_mesh, n=11):
         return as_mesh(fill_between(self, other_mesh=other_mesh, n=n))
 
     def flip(self, mask=None):
-        """Ensure positive cell volumes for `tria`, `tetra`, `quad` and `hexahedron`
+        r"""Ensure positive cell volumes for `tria`, `tetra`, `quad` and `hexahedron`
         cell types.
 
         Parameters
@@ -1080,8 +1080,8 @@ def flip(self, mask=None):
 
         Examples
         --------
-        A quad mesh with negative cell volumes occurs if one coordinate axis is multiplied
-        by -1. The error pops up if a region is created with this mesh.
+        A quad mesh with negative cell volumes occurs if one coordinate axis is
+        multiplied by -1. The error pops up if a region is created with this mesh.
 
         >>> import numpy as np
         >>> import felupe as fem
@@ -1090,8 +1090,8 @@ def flip(self, mask=None):
         >>> mesh.update(points=mesh.points * np.array([[-1, 1]]))
         >>> region = fem.RegionQuad(mesh)
 
-        The sum of the differential volumes :math:`V = \sum_c \sum_q dV_{qc}` is evaluated
-        to -1.0.
+        The sum of the differential volumes :math:`V = \sum_c \sum_q dV_{qc}` is
+        evaluated to -1.0.
 
         >>> region.dV.sum()
         -1.0

diff --git a/src/felupe/mesh/_tools.py b/src/felupe/mesh/_tools.py
@@ -696,7 +696,7 @@ def translate(points, cells, cell_type, move, axis):
 
 @mesh_or_data
 def flip(points, cells, cell_type, mask=None):
-    """Ensure positive cell volumes for `tria`, `tetra`, `quad` and `hexahedron` cell
+    r"""Ensure positive cell volumes for `tria`, `tetra`, `quad` and `hexahedron` cell
     types.
 
     Parameters

diff --git a/src/felupe/tools/_newton.py b/src/felupe/tools/_newton.py
@@ -98,6 +98,9 @@ def fun_items(items, x, parallel=False):
         # assemble vector
         r = body.assemble.vector(field=body.field, **kwargs)
 
+        if body.assemble.multiplier is not None:
+            r *= body.assemble.multiplier
+
         # check and reshape vector
         if r.shape != shape:
             r.resize(*shape)
@@ -122,6 +125,9 @@ def jac_items(items, x, parallel=False):
         # assemble matrix
         K = body.assemble.matrix(**kwargs)
 
+        if body.assemble.multiplier is not None:
+            K *= body.assemble.multiplier
+
         # check and reshape matrix
         if K.shape != matrix.shape:
             K.resize(*shape)