Repository: adtzlr/felupe

Issue:
Top-level Field for multiple solid bodies
Currently, it's not possible to combine multiple solid bodies with different element types. At least it is complicated.

# Tasks
- [x] Add a vertex element and a vertex-region. See #940
- [x] Add a method for a merged vertex mesh `MeshContainer.as_vertex_mesh()`. See #940 
- [x] Add a classmethod `Field.from_mesh_container(container)` or something similar. See #940
- [x] Handle mixed-fields (just add the additional fields to the top-level field) and take care of the offsets.

~The top-level field container should also handle different dimensions!?~ No.

Discussion:

Solution patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1bbf3676..20223e39 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,10 @@ All notable changes to this project will be documented in this file. The format
 - Add support for non-symmetric bilinear mixed forms in `IntegralForm`.
 - Add `element.Element` to the top-level package namespace.
 - Raise a TypeError if a mesh container is used as the `mesh`-argument in a region. The error message explains that the container is not supported. A mesh must be used instead.
+- Add the `Vertex` element formulation.
+- Add a vertex-region `RegionVertex`.
+- Add `MeshContainer.as_vertex_mesh()` to create a merged vertex mesh from the meshes of the mesh container.
+- Add `Field.from_mesh_container(container)` to create a top-level field on a vertex mesh.
 
 ### Changed
 - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`.
diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py
index a5e16468..e1b623a6 100644
--- a/src/felupe/__init__.py
+++ b/src/felupe/__init__.py
@@ -58,7 +58,7 @@ from .element import (
     TetraMINI,
     Triangle,
     TriangleMINI,
-    TriQuadraticHexahedron,
+    Vertex,
 )
 from .field import (
     Field,
@@ -113,6 +113,7 @@ from .region import (
     RegionTriangleMINI,
     RegionTriQuadraticHexahedron,
     RegionTriQuadraticHexahedronBoundary,
+    RegionVertex,
 )
 from .tools import hello_world, newtonrhapson, project, runs_on, save, topoints
 from .view import ViewField, ViewMesh
@@ -187,6 +188,7 @@ __all__ = [
     "Triangle",
     "TriangleMINI",
     "TriQuadraticHexahedron",
+    "Vertex",
     "Circle",
     "Cube",
     "Grid",
@@ -235,6 +237,7 @@ __all__ = [
     "RegionTriangleMINI",
     "RegionTriQuadraticHexahedron",
     "RegionTriQuadraticHexahedronBoundary",
+    "RegionVertex",
     "newtonrhapson",
     "project",
     "save",
diff --git a/src/felupe/element/__init__.py b/src/felupe/element/__init__.py
index 9ecbc8f7..461d54ff 100644
--- a/src/felupe/element/__init__.py
+++ b/src/felupe/element/__init__.py
@@ -15,6 +15,7 @@ from ._line import Line
 from ._quad import BiQuadraticQuad, ConstantQuad, Quad, QuadraticQuad
 from ._tetra import QuadraticTetra, Tetra, TetraMINI
 from ._triangle import QuadraticTriangle, Triangle, TriangleMINI
+from ._vertex import Vertex
 
 __all__ = [
     "Element",
@@ -37,4 +38,5 @@ __all__ = [
     "lagrange_line",
     "lagrange_quad",
     "lagrange_hexahedron",
+    "Vertex",
 ]
diff --git a/src/felupe/element/_vertex.py b/src/felupe/element/_vertex.py
new file mode 100644
index 00000000..b180c353
--- /dev/null
+++ b/src/felupe/element/_vertex.py
@@ -0,0 +1,47 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+from ._base import Element
+
+
+class Vertex(Element):
+    r"""A vertex element formulation with constant shape functions.
+
+    Notes
+    -----
+    The vertex element is defined by one point.
+    """
+
+    def __init__(self):
+        self.points = np.array([[0.0]])
+        self.cells = np.arange(len(self.points)).reshape(1, -1)
+        self.cell_type = "vertex"
+
+    def function(self, r):
+        "Return the shape functions at given coordinate (r)."
+        return np.ones(1)
+
+    def gradient(self, r):
+        "Return the gradient of shape functions at given coordinate (r)."
+        return np.zeros((1, 1))
+
+    def hessian(self, rs):
+        "Return the hessian of shape functions at given coordinate (r)."
+        return np.zeros((1, 1, 1))
diff --git a/src/felupe/field/_base.py b/src/felupe/field/_base.py
index 835f56b1..6b53714d 100644
--- a/src/felupe/field/_base.py
+++ b/src/felupe/field/_base.py
@@ -22,6 +22,7 @@ import numpy as np
 
 from ..math import identity
 from ..math import sym as symmetric
+from ..region import RegionVertex
 from ._container import FieldContainer
 from ._indices import Indices
 
@@ -138,6 +139,18 @@ class Field:
 
         return cai, ai
 
+    @classmethod
+    def from_mesh_container(cls, mesh_container, dim=None, values=0.0):
+        "Create a Field on a vertex mesh from a mesh container."
+
+        mesh = mesh_container.as_vertex_mesh()
+        region = RegionVertex(mesh)
+
+        if dim is None:
+            dim = mesh.dim
+
+        return cls(region, dim=dim, values=values)
+
     def grad(self, sym=False, dtype=None, out=None, order="C"):
         r"""Gradient as partial derivative of field values w.r.t. undeformed
         coordinates, evaluated at the integration points of all cells in the region.
diff --git a/src/felupe/mesh/_container.py b/src/felupe/mesh/_container.py
index 17baa2ad..8cf234ff 100644
--- a/src/felupe/mesh/_container.py
+++ b/src/felupe/mesh/_container.py
@@ -279,6 +279,11 @@ class MeshContainer:
 
         return meshio.Mesh(self.points, cells, **kwargs)
 
+    def as_vertex_mesh(self):
+        "Return a merged vertex-mesh."
+        cells = np.unique(np.concatenate([mesh.cells.ravel() for mesh in self.meshes]))
+        return Mesh(self.points, cells.reshape(-1, 1), cell_type="vertex")
+
     def copy(self):
         "Return a deepcopy of the mesh container."
         return deepcopy(self)
diff --git a/src/felupe/mesh/_convert.py b/src/felupe/mesh/_convert.py
index 4b978647..3c8377f7 100644
--- a/src/felupe/mesh/_convert.py
+++ b/src/felupe/mesh/_convert.py
@@ -59,6 +59,7 @@ def cell_types():
     import pyvista as pv
 
     cell_types = [
+        ("vertex", pv.CellType.VERTEX),
         ("line", pv.CellType.LINE),
         ("triangle", pv.CellType.TRIANGLE),
         ("triangle6", pv.CellType.QUADRATIC_TRIANGLE),
diff --git a/src/felupe/region/__init__.py b/src/felupe/region/__init__.py
index 8569957b..4086467b 100644
--- a/src/felupe/region/__init__.py
+++ b/src/felupe/region/__init__.py
@@ -22,6 +22,7 @@ from ._templates import (
     RegionTriangleMINI,
     RegionTriQuadraticHexahedron,
     RegionTriQuadraticHexahedronBoundary,
+    RegionVertex,
 )
 
 __all__ = [
@@ -48,4 +49,5 @@ __all__ = [
     "RegionTriangleMINI",
     "RegionTriQuadraticHexahedron",
     "RegionTriQuadraticHexahedronBoundary",
+    "RegionVertex",
 ]
diff --git a/src/felupe/region/_templates.py b/src/felupe/region/_templates.py
index 3cdac385..ec25f8cb 100644
--- a/src/felupe/region/_templates.py
+++ b/src/felupe/region/_templates.py
@@ -32,6 +32,7 @@ from ..element import (
     Triangle,
     TriangleMINI,
     TriQuadraticHexahedron,
+    Vertex,
 )
 from ..mesh import Mesh
 from ..quadrature import GaussLegendre, GaussLegendreBoundary
@@ -852,3 +853,13 @@ class RegionQuadraticTetra(Region):
     def __init__(self, mesh, quadrature=TetraQuadrature(order=2), grad=True, **kwargs):
         element = QuadraticTetra()
         super().__init__(mesh, element, quadrature, grad=grad, **kwargs)
+
+
+class RegionVertex(Region):
+    "A region with a vertex element."
+
+    def __init__(
+        self, mesh, quadrature=GaussLegendre(order=0, dim=1), grad=False, **kwargs
+    ):
+        element = Vertex()
+        super().__init__(mesh, element, quadrature, grad=grad, **kwargs)


Test patch:
diff --git a/tests/test_element.py b/tests/test_element.py
index a8c78571..afc92a55 100644
--- a/tests/test_element.py
+++ b/tests/test_element.py
@@ -30,6 +30,20 @@ import pytest
 import felupe as fem
 
 
+def test_vertex1():
+    vertex1 = fem.element.Vertex()
+
+    r = [0]
+
+    h = vertex1.function(r)
+    dhdr = vertex1.gradient(r)
+    d2hdrdr = vertex1.hessian(r)
+
+    assert h[0] == 1.0
+    assert np.all(dhdr[0] == 0.0)
+    assert np.all(d2hdrdr == 0.0)
+
+
 def test_line2():
     line2 = fem.element.Line()
 
@@ -185,6 +199,7 @@ def test_tri3():
     assert np.all(dhdr[0] == -1)
     assert np.all(d2hdrdr == 0)
 
+
 def test_tri6():
     tri6 = fem.element.QuadraticTriangle()
 
@@ -273,6 +288,8 @@ def test_aol():
 
 
 if __name__ == "__main__":
+    test_vertex1()
+
     test_line2()
     test_line_lagrange()
 
diff --git a/tests/test_field.py b/tests/test_field.py
index 7b6d0135..a3d8f2be 100644
--- a/tests/test_field.py
+++ b/tests/test_field.py
@@ -263,6 +263,34 @@ def test_link():
     assert field1[1].values is field2[1].values
 
 
+def test_toplevel():
+    meshes = [
+        fem.Rectangle(n=3),
+        fem.Rectangle(n=3).translate(1, axis=0).triangulate(),
+    ]
+    container = fem.MeshContainer(meshes, merge=True)
+    field = fem.Field.from_mesh_container(container).as_container()
+    regions = [
+        fem.RegionQuad(container.meshes[0]),
+        fem.RegionTriangle(container.meshes[1]),
+    ]
+    fields = [
+        fem.FieldContainer([fem.FieldPlaneStrain(regions[0], dim=2)]),
+        fem.FieldContainer([fem.FieldPlaneStrain(regions[1], dim=2)]),
+    ]
+    boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)
+    umats = [
+        fem.LinearElastic(E=2.1e5, nu=0.3),
+        fem.LinearElastic(E=1.0, nu=0.3),
+    ]
+    solids = [
+        fem.SolidBody(umat=umats[0], field=fields[0]),
+        fem.SolidBody(umat=umats[1], field=fields[1]),
+    ]
+    step = fem.Step(items=solids, boundaries=boundaries)
+    fem.Job(steps=[step]).evaluate(x0=field)
+
+
 if __name__ == "__main__":
     test_axi()
     test_3d()
@@ -270,3 +298,4 @@ if __name__ == "__main__":
     test_mixed_lagrange()
     test_view()
     test_link()
+    test_toplevel()
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 4c66604c..835cb286 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -374,6 +374,10 @@ def test_container():
     assert np.allclose(container_2[0].points, mesh_1.points)
     assert np.allclose(container_2[0].cells, mesh_1.cells)
 
+    mesh = container.as_vertex_mesh()
+
+    assert mesh.cells.shape[1] == 1
+
 
 def test_read(filename="tests/mesh.bdf"):
     mesh = fem.mesh.read(filename=filename, dim=2)[0]
diff --git a/tests/test_region.py b/tests/test_region.py
index b0cd7180..47bdb0a9 100644
--- a/tests/test_region.py
+++ b/tests/test_region.py
@@ -32,6 +32,12 @@ import felupe as fem
 
 
 def test_region():
+
+    mesh = fem.Point()
+    r = fem.RegionVertex(mesh)
+
+    assert r.h.flags["C_CONTIGUOUS"]
+
     mesh = fem.Rectangle()
     r = fem.RegionQuad(mesh, uniform=True)