adtzlr · adtzlr · Jun 17, 2024 · Jun 17, 2024 · Jun 17, 2024 · Jun 17, 2024
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -3,6 +3,10 @@ All notable changes to this project will be documented in this file. The format
 
 ## [Unreleased]
 
+### Added
+- Add `Region.astype(dtype=None)` to copy and cast the region arrays to a specified type.
+- Add `Field.extract(dtype=None)` to cast the extracted field gradient/interpolated values to a specified type.
+
 ## [8.8.0] - 2024-06-16
 
 ### Added

diff --git a/docs/howto/mixed.rst b/docs/howto/mixed.rst
@@ -33,7 +33,7 @@ Next, let's create a meshed cube for a Hood-Taylor element formulation. The fami
 
     displacement = fem.Field(region_q2,  dim=3)
     pressure     = fem.Field(region_p1, dim=1)
-    volumeratio  = fem.Field(region_p1, dim=1, values=1)
+    volumeratio  = fem.Field(region_p1, dim=1, values=1.0)
 
     field = fem.FieldContainer(fields=[displacement, pressure, volumeratio])
     solid = fem.SolidBody(umat=umat, field=field)

diff --git a/src/felupe/assembly/expression/_linear.py b/src/felupe/assembly/expression/_linear.py
@@ -62,7 +62,12 @@ def integrate(self, weakform, kwargs={}, parallel=False):
         """
 
         values = np.zeros(
-            (len(self.v.basis), self.v.basis.shape[-4], self.v.basis.shape[-2], self.v.field.region.mesh.ncells)
+            (
+                len(self.v.basis),
+                self.v.basis.shape[-4],
+                self.v.basis.shape[-2],
+                self.v.field.region.mesh.ncells,
+            )
         )
 
         if not parallel:

diff --git a/src/felupe/constitution/hyperelasticity/core/_neo_hooke_nearly_incompressible.py b/src/felupe/constitution/hyperelasticity/core/_neo_hooke_nearly_incompressible.py
@@ -323,7 +323,7 @@ def hessian(self, x, mu=None, bulk=None, out=None):
 
         A4 = out
         if A4 is None:
-            A4 = np.zeros((*F.shape[:2], *F.shape[:2], *F.shape[-2:]))
+            A4 = np.zeros((*F.shape[:2], *F.shape[:2], *F.shape[-2:]), dtype=F.dtype)
         else:
             np.multiply(A4, 0, out=A4)
 

diff --git a/src/felupe/field/_axi.py b/src/felupe/field/_axi.py
@@ -36,7 +36,9 @@ class FieldAxisymmetric(Field):
     values : float or array, optional
         A single value for all components of the field or an array of
         shape (region.mesh.npoints, dim)`. Default is 0.0.
-
+    dtype : data-type or None, optional
+        Data-type of the array containing the field values.
+
     Notes
     -----
     * component 1 =  axial component
@@ -75,18 +77,18 @@ class FieldAxisymmetric(Field):
 
     """
 
-    def __init__(self, region, dim=2, values=0.0):
+    def __init__(self, region, dim=2, values=0.0, dtype=None):
         # init base Field
-        super().__init__(region, dim=dim, values=values)
+        super().__init__(region, dim=dim, values=values, dtype=dtype)
 
         # create scalar-valued field of radial point values
-        self.scalar = Field(region, dim=1, values=region.mesh.points[:, 1])
+        self.scalar = Field(region, dim=1, values=region.mesh.points[:, 1], dtype=dtype)
 
         # interpolate radial point values to integration points of each cell
         # in the region
         self.radius = self.scalar.interpolate()
 
-    def _interpolate_2d(self, out=None):
+    def _interpolate_2d(self, dtype=None, out=None):
         """Interpolate 2D field values at points and evaluate them at the
         integration points of all cells in the region."""
 
@@ -97,18 +99,21 @@ def _interpolate_2d(self, out=None):
             "ca...,aqc->...qc",
             self.values[self.region.mesh.cells],
             self.region.h,
+            dtype=dtype,
             out=out,
         )
 
-    def interpolate(self, out=None):
+    def interpolate(self, dtype=None, out=None):
         # out-argument is not supported
         # if out is not None:
         #     out = out[:2]
 
         # extend dimension of in-plane 2d-gradient
-        return np.pad(self._interpolate_2d(out=None), ((0, 1), (0, 0), (0, 0)))
+        return np.pad(
+            self._interpolate_2d(dtype=dtype, out=None), ((0, 1), (0, 0), (0, 0))
+        )
 
-    def _grad_2d(self, sym=False, out=None):
+    def _grad_2d(self, sym=False, dtype=None, out=None):
         """In-plane 2D gradient as partial derivative of field values at points
         w.r.t. the undeformed coordinates, evaluated at the integration points
         of all cells in the region. Optionally, the symmetric part of the
@@ -118,6 +123,9 @@ def _grad_2d(self, sym=False, out=None):
         ---------
         sym : bool, optional (default is False)
             Calculate the symmetric part of the gradient.
+        dtype : data-type or None, optional
+            If provided, forces the calculation to use the data type specified. Default
+            is None.
         out : None or ndarray, optional
             A location into which the result is stored. If provided, it must have a
             shape that the inputs broadcast to. If not provided or None, a freshly-
@@ -138,6 +146,7 @@ def _grad_2d(self, sym=False, out=None):
             "ca...,aJqc->...Jqc",
             self.values[self.region.mesh.cells],
             self.region.dhdX,
+            dtype=dtype,
             out=out,
         )
 
@@ -146,7 +155,7 @@ def _grad_2d(self, sym=False, out=None):
         else:
             return g
 
-    def grad(self, sym=False, out=None):
+    def grad(self, sym=False, dtype=None, out=None):
         """3D-gradient as partial derivative of field values at points w.r.t.
         the undeformed coordinates, evaluated at the integration points of all
         cells in the region. Optionally, the symmetric part of the gradient is
@@ -162,6 +171,9 @@ def grad(self, sym=False, out=None):
         ---------
         sym : bool, optional
             Calculate the symmetric part of the gradient (default is False).
+        dtype : data-type or None, optional
+            If provided, forces the calculation to use the data type specified. Default
+            is None.
         out : None or ndarray, optional
             A location into which the result is stored. If provided, it must have a
             shape that the inputs broadcast to. If not provided or None, a freshly-
@@ -180,9 +192,12 @@ def grad(self, sym=False, out=None):
         #     out = out[:2, :2]
 
         # extend dimension of in-plane 2d-gradient
-        g = np.pad(self._grad_2d(sym=sym, out=None), ((0, 1), (0, 1), (0, 0), (0, 0)))
+        g = np.pad(
+            self._grad_2d(sym=sym, dtype=dtype, out=None),
+            ((0, 1), (0, 1), (0, 0), (0, 0)),
+        )
 
         # set dudX_33 = u_r / R
-        g[-1, -1] = self.interpolate()[1] / self.radius
+        g[-1, -1] = self.interpolate(dtype=dtype)[1] / self.radius
 
         return g
diff --git a/src/felupe/field/_base.py b/src/felupe/field/_base.py
@@ -39,6 +39,8 @@ class Field:
     values : float or array
         A single value for all components of the field or an array of shape
         `(region.mesh.npoints, dim)`. Default is 0.0.
+    dtype : data-type or None, optional
+        Data-type of the array containing the field values.
     **kwargs : dict, optional
         Extra class attributes for the field.
 
@@ -88,7 +90,7 @@ class Field:
 
     """
 
-    def __init__(self, region, dim=1, values=0.0, **kwargs):
+    def __init__(self, region, dim=1, values=0.0, dtype=None, **kwargs):
         self.region = region
         self.dim = dim
         self.shape = self.region.quadrature.npoints, self.region.mesh.ncells
@@ -107,8 +109,13 @@ def __init__(self, region, dim=1, values=0.0, **kwargs):
             else:
                 raise ValueError("Wrong shape of values.")
 
+            if dtype is not None:
+                self.values = self.values.astype(dtype)
+
         else:  # scalar value
-            self.values = np.ones((region.mesh.npoints, dim)) * values
+            self.values = np.full(
+                shape=(region.mesh.npoints, dim), fill_value=values, dtype=dtype
+            )
 
         eai, ai = self._indices_per_cell(self.region.mesh.cells, dim)
         self.indices = Indices(eai, ai, region, dim)
@@ -125,7 +132,7 @@ def _indices_per_cell(self, cells, dim):
 
         return cai, ai
 
-    def grad(self, sym=False, out=None):
+    def grad(self, sym=False, dtype=None, out=None):
         r"""Gradient as partial derivative of field values w.r.t. undeformed
         coordinates, evaluated at the integration points of all cells in the region.
         Optionally, the symmetric part the gradient is evaluated.
@@ -139,6 +146,9 @@ def grad(self, sym=False, out=None):
         ---------
         sym : bool, optional
             Calculate the symmetric part of the gradient (default is False).
+        dtype : data-type or None, optional
+            If provided, forces the calculation to use the data type specified. Default
+            is None.
         out : None or ndarray, optional
             A location into which the result is stored. If provided, it must have a
             shape that the inputs broadcast to. If not provided or None, a freshly-
@@ -159,15 +169,16 @@ def grad(self, sym=False, out=None):
             "ca...,aJqc->...Jqc",
             self.values[self.region.mesh.cells],
             self.region.dhdX,
+            dtype=dtype,
             out=out,
         )
 
         if sym:
-            return symmetric(g)
+            return symmetric(g, out=g)
         else:
             return g
 
-    def interpolate(self, out=None):
+    def interpolate(self, dtype=None, out=None):
         """Interpolate field values located at mesh-points to the quadrature points
         ``q`` of cells ``c`` in the region.
 
@@ -177,6 +188,9 @@ def interpolate(self, out=None):
 
         Arguments
         ---------
+        dtype : data-type or None, optional
+            If provided, forces the calculation to use the data type specified. Default
+            is None.
         out : None or ndarray, optional
             A location into which the result is stored. If provided, it must have a
             shape that the inputs broadcast to. If not provided or None, a freshly-
@@ -196,10 +210,11 @@ def interpolate(self, out=None):
             "ca...,aqc->...qc",
             self.values[self.region.mesh.cells],
             self.region.h,
+            dtype=dtype,
             out=None,
         )
 
-    def extract(self, grad=True, sym=False, add_identity=True, out=None):
+    def extract(self, grad=True, sym=False, add_identity=True, dtype=None, out=None):
         """Generalized extraction method which evaluates either the gradient or the
         field values at the integration points of all cells in the region. Optionally,
         the symmetric part of the gradient is evaluated and/or the identity matrix is
@@ -214,6 +229,9 @@ def extract(self, grad=True, sym=False, add_identity=True, out=None):
         add_identity : bool, optional
             Flag for the addition of the identity matrix
             if the gradient is evaluated (default is True).
+        dtype : data-type or None, optional
+            If provided, forces the calculation to use the data type specified. Default
+            is None.
         out : None or ndarray, optional
             A location into which the result is stored. If provided, it must have a
             shape that the inputs broadcast to. If not provided or None, a freshly-
@@ -232,17 +250,17 @@ def extract(self, grad=True, sym=False, add_identity=True, out=None):
         """
 
         if grad:
-            gr = self.grad(out=out)
+            gr = self.grad(out=out, dtype=dtype)
 
             if sym:
                 gr = symmetric(gr, out=gr)
 
             if add_identity:
-                gr = np.add(gr, identity(gr), out=gr)
+                gr = np.add(gr, identity(gr, dtype=dtype), out=gr)
 
             return gr
         else:
-            return self.interpolate(out=out)
+            return self.interpolate(out=out, dtype=dtype)
 
     def copy(self):
         "Return a copy of the field."

diff --git a/src/felupe/field/_container.py b/src/felupe/field/_container.py
@@ -111,7 +111,7 @@ def __repr__(self):
 
         return "\n".join([header, size, fields_header, *fields])
 
-    def extract(self, grad=True, sym=False, add_identity=True, out=None):
+    def extract(self, grad=True, sym=False, add_identity=True, dtype=None, out=None):
         """Generalized extraction method which evaluates either the gradient or the
         field values at the integration points of all cells in the region. Optionally,
         the symmetric part of the gradient is evaluated and/or the identity matrix is
@@ -129,6 +129,9 @@ def extract(self, grad=True, sym=False, add_identity=True, out=None):
         add_identity : bool, optional
             Flag for the addition of the identity matrix if the gradient is evaluated
             (default is True).
+        dtype : data-type or None, optional
+            If provided, forces the calculation to use the data type specified. Default
+            is None.
         out : None or ndarray, optional
             A location into which the result is stored. If provided, it must have a
             shape that the inputs broadcast to. If not provided or None, a freshly-
@@ -149,7 +152,7 @@ def extract(self, grad=True, sym=False, add_identity=True, out=None):
 
         grads = np.pad(grad, (0, len(self.fields) - 1))
         return tuple(
-            f.extract(g, sym, add_identity=add_identity, out=res)
+            f.extract(g, sym, add_identity=add_identity, dtype=dtype, out=res)
             for g, f, res in zip(grads, self.fields, out)
         )