Merge pull request #229 from adtzlr/unify-fields

Unify handling of `Field` and `FieldMixed`

README.md

@@ -26,17 +26,18 @@ region = fe.RegionHexahedron(fe.Cube(n=11))
 
 # add a displacement field and apply a uniaxial elongation on the cube
 displacement = fe.Field(region, dim=3)
-boundaries, dof0, dof1, ext0 = fe.dof.uniaxial(displacement, move=0.2, clamped=True)
+field = fe.FieldContainer([displacement])
+boundaries, dof0, dof1, offsets, ext0 = fe.dof.uniaxial(field, move=0.2, clamped=True)
 
 # define the constitutive material behavior and create a solid body
 umat = fe.NeoHooke(mu=1.0, bulk=2.0)
-solid = fe.SolidBody(umat, displacement)
+solid = fe.SolidBody(umat, field)
 
 # newton-rhapson procedure
-res = fe.newtonrhapson(bodies=[solid], dof1=dof1, dof0=dof0, ext0=ext0)
+res = fe.newtonrhapson(bodies=[solid], dof1=dof1, dof0=dof0, offsets=offsets, ext0=ext0)
 
 # save result
-fe.save(region, res.x, filename="result.vtk")
+fe.save(region, res.x, offsets=offsets, filename="result.vtk")
 ```
 
 <img src="https://raw.githubusercontent.com/adtzlr/felupe/main/docs/images/readme.png" width="600px"/>
@@ -53,6 +54,10 @@ All notable changes to this project will be documented in this file. The format
 - Add `SolidBody.evaluate.kirchhoff_stress()` method. Contrary to the Cauchy stress method, this gives correct results in incompressible plane stress.
 - Add `SolidBodyTensor` for tensor-based material definitions with state variables.
 - Add `bodies` argument to `newtonrhapson()`.
+- Add a container class for fields, `FieldContainer`.
+
+### Changed
+- Unify handling of `Field` and `FieldMixed`.
 
 ### Fixed
 - Fix `tovoigt()` helper for data with more or less than two trailing axes and 2D tensors.

docs/howto/axi.rst

@@ -10,18 +10,13 @@ For axisymmetric analyses an axisymmetric vector-valued field has to be created
     mesh = fe.Rectangle(n=3)
     region = fe.RegionQuad(mesh)
     u = fe.FieldAxisymmetric(region, dim=2)
+    field = fe.FieldContainer([u])
 
 Now it gets important: The 3x3 deformation gradient for an axisymmetric problem is obtained with :meth:`felupe.FieldAxisymmetric.grad` or :meth:`felupe.FieldAxisymmetric.extract` methods. For instances of :class:`felupe.FieldAxisymmetric` the gradient is modified to return a 3x3 gradient as described in :ref:`theory-axi`.
 
 ..  code-block:: python
 
-    dudX = fe.math.grad(u)
-
-or
-
-..  code-block:: python
-
-    F = u.extract(grad=True, sym=False, add_identity=True)
+    F = field.extract(grad=True, sym=False, add_identity=True)
 
 For simplicity, let's assume a (built-in) Neo-Hookean material.
 
@@ -30,13 +25,13 @@ For simplicity, let's assume a (built-in) Neo-Hookean material.
     umat = fe.NeoHooke(mu=1, bulk=5)
 
 
-FElupe provides an adopted :class:`felupe.IntegralFormAxisymmetric` class for the integration and the sparse matrix assemblage of axisymmetric problems. It uses the additional information (e.g. radial coordinates at integration points) stored in :class:`felupe.FieldAxisymmetric` to provide a consistent interface in comparison to default IntegralForms.
+FElupe provides an adopted :class:`felupe.IntegralFormAxisymmetric` class for the integration and the sparse matrix assemblage of axisymmetric problems under hood. It uses the additional information (e.g. radial coordinates at integration points) stored in :class:`felupe.FieldAxisymmetric` to provide a consistent interface in comparison to default IntegralForms.
 
 ..  code-block:: python
 
-    F = fe.math.identity(dudX) + dudX
+    dA = region.dV
 
-    r = fe.IntegralFormAxisymmetric(umat.gradient(F), u, dA).assemble()
-    K = fe.IntegralFormAxisymmetric(umat.hessian(F), u, dA).assemble()
+    r = fe.IntegralForm(umat.gradient(F), field, dA).assemble()
+    K = fe.IntegralForm(umat.hessian(F), field, dA, field).assemble()
 
-To sum up, for axisymmetric problems use :class:`felupe.FieldAxisymmetric` in conjunction with :class:`felupe.IntegralFormAxisymmetric`. Of course, Mixed-field formulations may be applied on axisymmetric scenarios too.
+To sum up, for axisymmetric problems use :class:`felupe.FieldAxisymmetric`. Of course, Mixed-field formulations may be applied on axisymmetric scenarios too.

docs/howto/composite.rst

@@ -12,6 +12,7 @@ This section demonstrates how to set up a problem with two regions, each associa
     mesh = fe.Cube(n=n)
     region = fe.RegionHexahedron(mesh)
     displacement = fe.Field(region, dim=3)
+    field = fe.FieldContainer([displacement])
 
 
 In a second step, sub-sets for points and cells are created from which two sub-regions and sub-fields are initiated. All field values are linked, that means they share their values array.
@@ -36,9 +37,11 @@ In a second step, sub-sets for points and cells are created from which two sub-r
     
     region_rubber = fe.RegionHexahedron(mesh_rubber)
     displacement_rubber = fe.Field(region_rubber, dim=3)
+    field_rubber = fe.FieldContainer([displacement_rubber])
 
     region_steel = fe.RegionHexahedron(mesh_steel)
     displacement_steel = fe.Field(region_steel, dim=3)
+    field_steel = fe.FieldContainer([displacement_steel])
 
     # link fields
     displacement_steel.values = displacement_rubber.values = displacement.values
@@ -49,7 +52,7 @@ The displacement boundaries are created on the total field.
 ..  code-block:: python
 
     boundaries, dof0, dof1, ext0 = fe.dof.uniaxial(
-        displacement, move=-0.25
+        field, move=-0.25
     )
 
 
@@ -60,8 +63,8 @@ The rubber is associated to a Neo-Hookean material formulation whereas the steel
     neohooke = fe.NeoHooke(mu=1.0, bulk=2.0)
     linearelastic = fe.LinearElastic(E=210000.0, nu=0.3)
 
-    rubber = fe.SolidBody(neohooke, displacement_rubber)
-    steel = fe.SolidBody(linearelastic, displacement_steel)
+    rubber = fe.SolidBody(neohooke, field_rubber)
+    steel = fe.SolidBody(linearelastic, field_steel)
 
 
 Inside the Newton-Rhapson iterations both the internal force vector and the tangent stiffness matrix are assembled and summed up from contributions of both solid bodies.
@@ -76,15 +79,15 @@ Inside the Newton-Rhapson iterations both the internal force vector and the tang
         K = rubber.assemble.matrix()
         K+= steel.assemble.matrix()
 
-        system = fe.solve.partition(displacement, K, dof1, dof0, r)
-        du = fe.solve.solve(*system, ext0)
+        system = fe.solve.partition(field, K, dof1, dof0, r)
+        dfield = np.split(fe.solve.solve(*system, ext0), field.offsets)
 
-        displacement += du
+        field += dfield
         
-        r = rubber.assemble.vector(displacement_rubber)
-        r+= steel.assemble.vector(displacement_steel)
+        r = rubber.assemble.vector(field_rubber)
+        r+= steel.assemble.vector(field_steel)
 
-        norm = fe.math.norm(du)
+        norm = fe.math.norm(dfield[0])
         print(iteration, norm)
 
         if norm < 1e-12:
@@ -109,10 +112,12 @@ Results and may be exported either for the total region or with stresses for sub
     s = rubber.evaluate.cauchy_stress()
     cauchy_stress = fe.project(fe.math.tovoigt(s), region_rubber)
     
-    fe.save(region, displacement, filename="result.vtk")
+    fe.save(region, field, filename="result.vtk")
 
-    fe.save(region_rubber, displacement_rubber, filename="result_rubber.vtk",
-        point_data={"CauchyStress": cauchy_stress})
+    fe.save(region_rubber, field_rubber,
+        filename="result_rubber.vtk", 
+        point_data={"CauchyStress": cauchy_stress}
+    )
 
 .. image:: images/composite_rubber_cauchy.png
    :width: 600px

docs/howto/forms.rst

@@ -1,7 +1,7 @@
 Linear and Bilinear Forms
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 
-FElupe requires a pre-evaluated array for the definition of a bilinear :class:`felupe.IntegralForm` object on interpolated field values or their gradients. While this has two benefits, namely a fast integration of the form is easy to code and the array may be computed in any programming language, sometimes numeric representations of analytic linear and bilinear form expressions may be easier in user-code and less error prone compared to the calculation of explicit second or fourth-order tensors. Therefore, FElupe provides a function decorator :func:`felupe.Form` as an easy-to-use high-level interface, similar to what `scikit-fem <https://github.com/kinnala/scikit-fem>`_ offers. While the :func:`felupe.Form` decorator handles both single and mixed fields, additional access to the underlying form objects is enabled by :class:`felupe.LinearForm` and :class:`felupe.BilinearForm` for single-field as well as :class:`felupe.LinearFormMixed` and :class:`felupe.BilinearFormMixed` for mixed-field problems. All these linear and bilinear form classes are similar, but not identical in their usage compared to :class:`felupe.IntegralForm`. They require a callable function (with optional arguments and keyword arguments) instead of a pre-computed array to be passed. The bilinear form of linear elasticity serves as a reference example for the demonstration on how to use this feature of FElupe. The stiffness matrix is assembled for a unit cube out of hexahedrons.
+FElupe requires a pre-evaluated array for the definition of a bilinear :class:`felupe.IntegralForm` object on interpolated field values or their gradients. While this has two benefits, namely a fast integration of the form is easy to code and the array may be computed in any programming language, sometimes numeric representations of analytic linear and bilinear form expressions may be easier in user-code and less error prone compared to the calculation of explicit second or fourth-order tensors. Therefore, FElupe provides a function decorator :func:`felupe.Form` as an easy-to-use high-level interface, similar to what `scikit-fem <https://github.com/kinnala/scikit-fem>`_ offers. The :func:`felupe.Form` decorator handles a field container. The form class is similar, but not identical in its usage compared to :class:`felupe.IntegralForm`. It requires a callable function (with optional arguments and keyword arguments) instead of a pre-computed array to be passed. The bilinear form of linear elasticity serves as a reference example for the demonstration on how to use this feature of FElupe. The stiffness matrix is assembled for a unit cube out of hexahedrons.
 
 ..  code-block:: python
 
@@ -10,6 +10,7 @@ FElupe requires a pre-evaluated array for the definition of a bilinear :class:`f
     mesh = fe.Cube(n=11)
     region = fe.RegionHexahedron(mesh)
     displacement = fe.Field(region, dim=3)
+    field = fe.FieldContainer([displacement])
 
 The bilinear form of linear elasticity is defined as
 
@@ -31,11 +32,16 @@ and implemented in FElupe closely to the analytic expression. The first two argu
 
     from felupe.math import ddot, trace, sym
     
-    @fe.Form(v=displacement, u=displacement, grad_v=True, grad_u=True, kwargs={"mu": 1.0, "lmbda": 2.0})
-    def linear_elasticity(gradv, gradu, mu, lmbda):
-        "Linear elasticity."
+    @fe.Form(v=field, u=field, grad_v=[True], grad_u=[True], kwargs={"mu": 1.0, "lmbda": 2.0})
+    def bilinearform():
+        "A container for a bilinear form."
         
-        de, e = sym(gradv), sym(gradu)
-        return 2 * mu * ddot(de, e) + lmbda * trace(de) * trace(e)
+        def linear_elasticity(gradv, gradu, mu, lmbda):
+            "Linear elasticity."
+        
+            de, e = sym(gradv), sym(gradu)
+            return 2 * mu * ddot(de, e) + lmbda * trace(de) * trace(e)
+        
+        return [linear_elasticity,]
 
-    K = linear_elasticity.assemble(v=displacement, u=displacement, parallel=False)
+    K = bilinearform.assemble(v=field, u=field, parallel=False)

docs/howto/mixed.rst

@@ -25,7 +25,7 @@ Next, let's create a meshed cube. Two regions, one for the displacements and ano
     pressure     = fe.Field(region0, dim=1)
     volumeratio  = fe.Field(region0, dim=1, values=1)
 
-    fields = fe.FieldMixed((displacement, pressure, volumeratio))
+    field = fe.FieldContainer((displacement, pressure, volumeratio))
 
 Boundary conditions are enforced in the same way as in Getting Started.
 
@@ -39,35 +39,35 @@ Boundary conditions are enforced in the same way as in Getting Started.
     boundaries["right"] = fe.Boundary(displacement, fx=f1, skip=(1, 0, 0))
     boundaries["move" ] = fe.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=-0.4)
 
-    dof0, dof1, offsets = fe.dof.partition(fields, boundaries)
-    ext0 = fe.dof.apply(fields, boundaries, dof0, offsets)
+    dof0, dof1 = fe.dof.partition(field, boundaries)
+    ext0 = fe.dof.apply(field, boundaries, dof0)
 
 The Newton-Rhapson iterations are coded quite similar to the one used in :ref:`tutorial-getting-started`. FElupe provides a Mixed-field version of it's :class:`felupe.IntegralForm`, called :class:`felupe.IntegralFormMixed`. It assumes that the first field operates on the gradient and all the others don't. The resulting system vector with incremental values of the fields has to be splitted at the field-offsets in order to update the fields.
 
 ..  code-block:: python
 
     for iteration in range(8):
 
-        F, p, J = fields.extract()
+        F, p, J = field.extract()
         
-        linearform   = fe.IntegralFormMixed(umat.gradient(F, p, J), fields, dV)
-        bilinearform = fe.IntegralFormMixed(umat.hessian(F, p, J), fields, dV, fields)
+        linearform   = fe.IntegralForm(umat.gradient([F, p, J]), field, dV)
+        bilinearform = fe.IntegralForm(umat.hessian([F, p, J]), field, dV, field)
 
         r = linearform.assemble().toarray()[:, 0]
         K = bilinearform.assemble()
         
-        system = fe.solve.partition(fields, K, dof1, dof0, r)
-        dfields = np.split(fe.solve.solve(*system, ext0), offsets)
+        system = fe.solve.partition(field, K, dof1, dof0, r)
+        dfield = np.split(fe.solve.solve(*system, ext0), field.offsets)
         
-        fields += dfields
+        field += dfield
 
-        norm = np.linalg.norm(dfields[0])
+        norm = np.linalg.norm(dfield[0])
         print(iteration, norm)
 
         if norm < 1e-12:
             break
 
-    fe.tools.save(region, fields, offsets=offsets, filename="result.vtk")
+    fe.tools.save(region, field, filename="result.vtk")
 
 The deformed cube is finally visualized by a VTK output file with the help of Paraview.
 

docs/howto/rbe2.rst

@@ -9,8 +9,8 @@ This How-To demonstrates the usage of multi-point constraints (also called MPC o
     import felupe as fe
 
     # mesh with one additional rbe2-control point
-    mesh = fe.Rectangle(n=11)
-    mesh.points = np.vstack((mesh.points, [2.0, 0.0]))
+    mesh = fe.Cube(n=11)
+    mesh.points = np.vstack((mesh.points, [2.0, 0.0, 0.0]))
     mesh.update(mesh.cells)
 
 An instance of :class:`felupe.MultiPointConstraint` defines the multi-point constraint. This instance provides two methods, :meth:`felupe.MultiPointConstraint.stiffness` and :meth:`felupe.MultiPointConstraint.residuals`.
@@ -28,12 +28,15 @@ Finally, add the results of these methods to the internal force vector or the st
 
 ..  code-block:: python
 
-    # ...
+    region = felupe.RegionHexahedron(mesh)
+    displacement = felupe.Field(region, dim=3)
+    field = felupe.FieldContainer([displacement])
+    umat = felupe.constitution.NeoHooke(mu=1.0, bulk=2.0)
 
     K = fe.IntegralForm(
-        umat.hessian(field.extract()), field, region.dV, field, True, True
+        umat.hessian(field.extract()), field, region.dV, field
     ).assemble() + MPC.stiffness()
 
     r = fe.IntegralForm(
-        umat.gradient(field.extract()), field, region.dV, None, True
+        umat.gradient(field.extract()), field, region.dV
     ).assemble() + MPC.residuals(field)

docs/howto/solid.rst

@@ -16,17 +16,18 @@ The generation of internal force vectors or stiffness matrices of solid bodies a
     mesh = fe.Cube(n=6)
     region = fe.RegionHexahedron(mesh)
     displacement = fe.Field(region, dim=3)
+    field = fe.FieldContainer([displacement])
     
-    body = fe.SolidBody(umat=neohooke, field=displacement)
-    internal_force = body.assemble.vector(displacement, parallel=False, jit=False)
-    stiffness_matrix = body.assemble.matrix(displacement, parallel=False, jit=False)
+    body = fe.SolidBody(umat=neohooke, field=field)
+    internal_force = body.assemble.vector(field, parallel=False, jit=False)
+    stiffness_matrix = body.assemble.matrix(field, parallel=False, jit=False)
 
 
 During assembly, several results are stored, e.g. the gradient of the strain energy density function per unit undeformed volume (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.
 
 ..  code-block:: python
     
-    F = body.results.kinematics[0]
+    F = body.results.kinematics
     P = body.results.stress
     A = body.results.elasticity
 
@@ -35,9 +36,9 @@ The Cauchy stress tensor, as well as the gradient and the hessian of the strain
 
 ..  code-block:: python
     
-    P = body.evaluate.gradient(displacement)
-    A = body.evaluate.hessian(displacement)
-    s = body.evaluate.cauchy_stress(displacement)
+    P = body.evaluate.gradient(field)
+    A = body.evaluate.hessian(field)
+    s = body.evaluate.cauchy_stress(field)
 
 
 Pressure Boundary on a Solid Body
@@ -53,17 +54,18 @@ The generation of internal force vectors or stiffness matrices of pressure bound
         mask=mesh.points[:, 0] == 0, # select a subset of faces on the surface
     )
     
-    displacement_boundary = fe.Field(region_pressure)
+    displacement_boundary = fe.Field(region_pressure, dim=3)
+    field_boundary = fe.FieldContainer([displacement_boundary])
     displacement_boundary.values = displacement.values # link field values
     
-    body_pressure = fe.SolidBodyPressure(field=displacement)
+    body_pressure = fe.SolidBodyPressure(field=field_boundary)
     
-    internal_force_pressure = body.assemble.vector(
-        field=displacement, parallel=False, jit=False, resize=internal_force
+    internal_force_pressure = body_pressure.assemble.vector(
+        field=field_boundary, parallel=False, jit=False, resize=internal_force
     )
     
-    stiffness_matrix_pressure = body.assemble.matrix(
-        field=displacement, parallel=False, jit=False, resize=stiffness_matrix
+    stiffness_matrix_pressure = body_pressure.assemble.matrix(
+        field=field_boundary, parallel=False, jit=False, resize=stiffness_matrix
     )
 
 
@@ -78,7 +80,7 @@ For axisymmetric problems the boundary region has to be created with the attribu
         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,
+        ensure_3d=True, # flag for axisymmetric boundary region
     )
     
     displacement = fe.FieldAxisymmetric(region)

docs/tutorial.rst

@@ -12,4 +12,3 @@ This section is all about learning. Each tutorial focusses on a some lessons to
    tutorial/platewithhole
    tutorial/beam
    tutorial/shear
-   tutorial/numcont