diff --git a/CHANGELOG.md b/CHANGELOG.md
index e04dafda..90b8d656 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,10 +14,12 @@ All notable changes to this project will be documented in this file. The format
 - Add `constitution.jax.Material(..., jacobian=None)` with JAX as backend. A custom jacobian-callable may be passed to switch between forward- and backward-mode automatic differentiation.
 - Add material models for JAX-based materials: `felupe.constitution.jax.models.hyperelastic.mooney_rivlin()`, `felupe.constitution.jax.models.hyperelastic.yeoh()`, `felupe.constitution.jax.models.hyperelastic.third_order_deformation()`, `felupe.constitution.jax.models.lagrange.morph()`, `felupe.constitution.jax.models.lagrange.morph_representative_directions()`.
 - Add `felupe.constitution.jax.total_lagrange()`, `felupe.constitution.jax.updated_lagrange()` and `felupe.constitution.jax.isochoric_volumetric_split()` function decorators for the JAX hyperelastic material class.
+- Add an optional keyword-argument `Region.astype(copy=True)` to modify the data types of the arrays of the region in-place if `copy=False`.
 
 ### Changed
 - Change default `np.einsum(..., order="K")` to `np.einsum(..., order="C")` in the methods of `Field`, `FieldAxisymmetric`, `FieldPlaneStrain` and `FieldContainer`.
 - Change supported Python versions to 3.9 - 3.12.
+- Change the `dtype`-argument in `Region.astype(dtype)` from an optional to a required argument.
 
 ### Fixed
 - Fix the number of points for non-disconnected dual meshes. This reduces the assembled (sparse) vector- and matrix-shapes, which are defined on mixed-fields.
diff --git a/src/felupe/region/_region.py b/src/felupe/region/_region.py
index cc088cc5..dd12b1ea 100644
--- a/src/felupe/region/_region.py
+++ b/src/felupe/region/_region.py
@@ -149,14 +149,17 @@ def __init__(self, mesh, element, quadrature, grad=True, hess=False, uniform=Fal
         self.evaluate_hessian = hess
         self.reload(mesh=mesh, element=element, quadrature=quadrature, uniform=uniform)
 
-    def astype(self, dtype=None):
+    def astype(self, dtype, copy=True):
         """Copy the region and cast the arrays to a specified type.
 
         Parameters
         ----------
-        dtype : data-type or None, optional
-            The data-type of the arrays of the Region. If None, a copy of the Region is
-            returned.
+        dtype : str or dtype
+            Typecode or data-type to which the arrays of the region are cast.
+        copy : bool, optional
+            By default, astype always returns a copy of the region with newly allocated
+            arrays. If False, the arrays of the input region are modified and the input
+            region is returned. Default is True.
 
         Returns
         -------
@@ -168,7 +171,11 @@ def astype(self, dtype=None):
         felupe.region.copy : Return a copy of the region and reload it if necessary.
         """
 
-        region = self.copy(uniform=self.uniform)
+        region = self
+
+        if copy:
+            region = region.copy(uniform=self.uniform)
+
         region.h = region.h.astype(dtype)
         region.dhdr = region.dhdr.astype(dtype)
 
@@ -178,6 +185,10 @@ def astype(self, dtype=None):
             region.dhdX = region.dhdX.astype(dtype)
             region.dV = region.dV.astype(dtype)
 
+        if region.evaluate_hessian:
+            region.d2hdrdr = region.d2hdrdr.astype(dtype)
+            region.d2hdXdX = region.d2hdXdX.astype(dtype)
+
         return region
 
     def copy(
diff --git a/tests/test_dtype.py b/tests/test_dtype.py
index 273c2126..08ff2862 100644
--- a/tests/test_dtype.py
+++ b/tests/test_dtype.py
@@ -32,7 +32,7 @@
 
 def test_dtype(dtype=np.float32, tol=1e-3):
     mesh = fem.Cube(n=3)
-    region = fem.RegionHexahedron(mesh).astype(dtype)
+    region = fem.RegionHexahedron(mesh).astype(dtype, copy=False)
     displacement = fem.Field(region, dim=3, dtype=dtype)
     field = fem.FieldContainer([displacement])
 
@@ -62,7 +62,7 @@ def test_dtype_axi():
 
 def test_dtype_planestrain():
     mesh = fem.Rectangle(n=3)
-    region = fem.RegionQuad(mesh).astype(np.float32)
+    region = fem.RegionQuad(mesh, hess=True).astype(np.float32)
     displacement = fem.FieldPlaneStrain(region, dtype=np.float32)
     field = fem.FieldContainer([displacement])