diff --git a/src/Thunderbolt.jl b/src/Thunderbolt.jl
index 780446dea..fd6b20e11 100644
--- a/src/Thunderbolt.jl
+++ b/src/Thunderbolt.jl
@@ -15,52 +15,38 @@ import Ferrite: vertices, edges, faces, sortedge, sortface
 
 import Krylov: CgSolver
 
-include("collections.jl")
 include("utils.jl")
 
 include("mesh/meshes.jl")
-include("mesh/coordinate_systems.jl")
-include("mesh/tools.jl")
-include("mesh/generators.jl")
 
-
-# TODO where to put these?
-abstract type SteadyStateInternalVariable end
-
-include("modeling/coefficients.jl")
-
-include("modeling/boundary_conditions.jl")
+# Note that some modules below have an "interface.jl" but this one has only a "common.jl".
+# This is simply because there is no modeling interface, but just individual physics modules and couplers.
+include("modeling/common.jl")
 
 include("modeling/microstructure.jl")
 
 include("modeling/electrophysiology.jl")
+include("modeling/solid_mechanics.jl")
+include("modeling/fluid_mechanics.jl")
 
-include("modeling/solid/energies.jl")
-include("modeling/solid/contraction.jl")
-include("modeling/solid/active.jl")
-include("modeling/solid/drivers.jl") # TODO better name. This is basically the quadrature point routine.
+include("modeling/multiphysics.jl")
 
-include("modeling/fluid/lumped.jl")
-
-include("modeling/coupler/interface.jl")
-include("modeling/coupler/fsi.jl")
-
-include("modeling/problems.jl")
+include("modeling/problems.jl") # This is not really "modeling" but a glue layer to translate from model to solver via a discretization
 
 include("solver/interface.jl")
 include("solver/operator.jl")
 include("solver/newton_raphson.jl")
 include("solver/load_stepping.jl")
-include("solver/backward_euler.jl")
+include("solver/euler.jl")
 include("solver/partitioned_solver.jl")
 include("solver/operator_splitting.jl")
 
 include("discretization/interface.jl")
 include("discretization/fem.jl")
 
-
 include("io.jl")
 
+# TODO put exports into the individual submodules above!
 export
     solve,
     # Coefficients
diff --git a/src/mesh/meshes.jl b/src/mesh/meshes.jl
index 6fde629ae..e7046556a 100644
--- a/src/mesh/meshes.jl
+++ b/src/mesh/meshes.jl
@@ -11,143 +11,13 @@
 #     defining_nodes::NTuple{3,Int}
 # end
 
+# TODO we might want to add this to Ferrite (and especially FerriteViz) in one or another way. Maybe traits are better, because they allow more extensibility.
 const LinearCellGeometry = Union{Hexahedron, Tetrahedron, Pyramid, Wedge, Triangle, Quadrilateral, Line}
 
 elementtypes(grid::Grid{3,Hexahedron}) = @SVector [Hexahedron]
 elementtypes(grid::Grid{3,Tetrahedron}) = @SVector [Tetrahedron]
 
-"""
-    SimpleMesh3D{C <: AbstractCell, T <: Real} <: AbstractGrid{3}
-
-A grid which also has information abouts its vertices, faces and edges.
-"""
-struct SimpleMesh3D{C <: AbstractCell, T <: Real} <: AbstractGrid{3}
-    grid::Grid{3, C, T}
-    mfaces::OrderedDict{NTuple{3,Int}, Int} # Maps "sortface"-representation to id
-    medges::OrderedDict{NTuple{2,Int}, Int} # Maps "sortedge"-representation to id
-    mvertices::OrderedDict{Int, Int} # Maps node to id
-    number_of_cells_by_type::Dict{DataType, Int}
-end
-
-global_edges(mgrid::SimpleMesh3D, cell) = [mgrid.medges[sedge] for sedge ∈ first.(sortedge.(edges(cell)))]
-global_faces(mgrid::SimpleMesh3D, cell) = [mgrid.mfaces[sface] for sface ∈ first.(sortface.(faces(cell)))]
-global_vertices(mgrid::SimpleMesh3D, cell) = [mgrid.mvertices[v] for v ∈ vertices(cell)]
-
-num_nodes(mgrid::SimpleMesh3D) = length(mgrid.grid.nodes)
-num_faces(mgrid::SimpleMesh3D) = length(mgrid.mfaces)
-num_edges(mgrid::SimpleMesh3D) = length(mgrid.medges)
-num_vertices(mgrid::SimpleMesh3D) = length(mgrid.mvertices)
-
-elementtypes(::SimpleMesh3D{Tetrahedron}) = @SVector [Tetrahedron]
-elementtypes(::SimpleMesh3D{Hexahedron}) = @SVector [Hexahedron]
-
-function to_mesh(grid::Grid{3,C,T}) where {C, T}
-    mfaces = OrderedDict{NTuple{3,Int}, Int}()
-    medges = OrderedDict{NTuple{2,Int}, Int}()
-    mvertices = OrderedDict{Int, Int}()
-    next_face_idx = 1
-    next_edge_idx = 1
-    next_vertex_idx = 1
-    number_of_cells_by_type = Dict{DataType, Int}()
-    for cell ∈ getcells(grid)
-        cell_type = typeof(cell)
-        if haskey(number_of_cells_by_type, cell_type)
-            number_of_cells_by_type[cell_type] += 1
-        else
-            number_of_cells_by_type[cell_type] = 1
-        end
-
-        for v ∈ vertices(cell)
-            if !haskey(mvertices, v)
-                mvertices[v] = next_vertex_idx
-                next_vertex_idx += 1
-            end
-        end
-        for e ∈ first.(sortedge.(edges(cell)))
-            if !haskey(medges, e)
-                medges[e] = next_edge_idx
-                next_edge_idx += 1
-            end
-        end
-        for f ∈ first.(sortface.(faces(cell)))
-            if !haskey(mfaces, f)
-                mfaces[f] = next_face_idx
-                next_face_idx += 1
-            end
-        end
-    end
-    return SimpleMesh3D(grid, mfaces, medges, mvertices, number_of_cells_by_type)
-end
-
-
-"""
-    SimpleMesh2D{C <: AbstractCell, T <: Real} <: AbstractGrid{2}
-
-A grid which also has information abouts its vertices, faces and edges.
-"""
-struct SimpleMesh2D{C <: AbstractCell, T <: Real} <: AbstractGrid{2}
-    grid::Grid{2, C, T}
-    mfaces::OrderedDict{NTuple{2,Int}, Int} # Maps "sortface"-representation to id
-    mvertices::OrderedDict{Int, Int} # Maps node to id
-    number_of_cells_by_type::Dict{DataType, Int}
-end
-
-function to_mesh(grid::Grid{2,C,T}) where {C, T}
-    mfaces = OrderedDict{NTuple{2,Int}, Int}()
-    mvertices = OrderedDict{Int, Int}()
-    next_face_idx = 1
-    next_vertex_idx = 1
-    number_of_cells_by_type = Dict{DataType, Int}()
-    for cell ∈ getcells(grid)
-        cell_type = typeof(cell)
-        if haskey(number_of_cells_by_type, cell_type)
-            number_of_cells_by_type[cell_type] += 1
-        else
-            number_of_cells_by_type[cell_type] = 1
-        end
-
-        for v ∈ vertices(cell)
-            if !haskey(mvertices, v)
-                mvertices[v] = next_vertex_idx
-                next_vertex_idx += 1
-            end
-        end
-        for f ∈ first.(sortface.(faces(cell)))
-            if !haskey(mfaces, f)
-                mfaces[f] = next_face_idx
-                next_face_idx += 1
-            end
-        end
-    end
-    return SimpleMesh2D(grid, mfaces, mvertices, number_of_cells_by_type)
-end
-
-
-global_faces(mgrid::SimpleMesh2D, cell) = [mgrid.mfaces[sface] for sface ∈ first.(sortface.(faces(cell)))]
-global_vertices(mgrid::SimpleMesh2D, cell) = [mgrid.mvertices[v] for v ∈ vertices(cell)]
-
-num_nodes(mgrid::SimpleMesh2D) = length(mgrid.grid.nodes)
-num_faces(mgrid::SimpleMesh2D) = length(mgrid.mfaces)
-num_vertices(mgrid::SimpleMesh2D) = length(mgrid.mvertices)
-
-elementtypes(::SimpleMesh2D{Triangle}) = @SVector [Triangle]
-elementtypes(::SimpleMesh2D{Quadrilateral}) = @SVector [Quadrilateral]
-
-# Ferrite compat layer for the mesh
-@inline Ferrite.getncells(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getncells(mesh.grid)
-@inline Ferrite.getcells(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getcells(mesh.grid)
-@inline Ferrite.getcells(mesh::Union{SimpleMesh2D,SimpleMesh3D}, v::Union{Int, Vector{Int}}) = Ferrite.getcells(mesh.grid, v)
-@inline Ferrite.getcells(mesh::Union{SimpleMesh2D,SimpleMesh3D}, setname::String) = Ferrite.getcells(mesh.grid, setname)
-@inline Ferrite.getcelltype(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getcelltype(mesh.grid)
-@inline Ferrite.getcelltype(mesh::Union{SimpleMesh2D,SimpleMesh3D}, i::Int) = Ferrite.getcelltype(mesh.grid)
-
-@inline Ferrite.getnnodes(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getnnodes(mesh.grid)
-@inline Ferrite.nnodes_per_cell(mesh::Union{SimpleMesh2D,SimpleMesh3D}, i::Int=1) = Ferrite.nnodes_per_cell(mesh.grid, i)
-@inline Ferrite.getnodes(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getnodes(mesh.grid)
-@inline Ferrite.getnodes(mesh::Union{SimpleMesh2D,SimpleMesh3D}, v::Union{Int, Vector{Int}}) = Ferrite.getnodes(mesh.grid, v)
-@inline Ferrite.getnodes(mesh::Union{SimpleMesh2D,SimpleMesh3D}, setname::String) = Ferrite.getnodes(mesh.grid, setname)
-
-@inline Ferrite.vtk_grid(filename::AbstractString, mesh::Union{SimpleMesh2D,SimpleMesh3D}; kwargs...) = Ferrite.vtk_grid(filename, mesh.grid, kwargs...)
-
-@inline Ferrite.get_coordinate_type(::SimpleMesh2D{C,T}) where {C,T} = Vec{2,T} 
-@inline Ferrite.get_coordinate_type(::SimpleMesh3D{C,T}) where {C,T} = Vec{3,T} 
+include("simple_meshes.jl")
+include("coordinate_systems.jl")
+include("tools.jl")
+include("generators.jl")
diff --git a/src/mesh/simple_meshes.jl b/src/mesh/simple_meshes.jl
new file mode 100644
index 000000000..35574d98e
--- /dev/null
+++ b/src/mesh/simple_meshes.jl
@@ -0,0 +1,136 @@
+
+"""
+SimpleMesh3D{C <: AbstractCell, T <: Real} <: AbstractGrid{3}
+
+A grid which also has information abouts its vertices, faces and edges.
+"""
+struct SimpleMesh3D{C <: AbstractCell, T <: Real} <: AbstractGrid{3}
+grid::Grid{3, C, T}
+mfaces::OrderedDict{NTuple{3,Int}, Int} # Maps "sortface"-representation to id
+medges::OrderedDict{NTuple{2,Int}, Int} # Maps "sortedge"-representation to id
+mvertices::OrderedDict{Int, Int} # Maps node to id
+number_of_cells_by_type::Dict{DataType, Int}
+end
+
+global_edges(mgrid::SimpleMesh3D, cell) = [mgrid.medges[sedge] for sedge ∈ first.(sortedge.(edges(cell)))]
+global_faces(mgrid::SimpleMesh3D, cell) = [mgrid.mfaces[sface] for sface ∈ first.(sortface.(faces(cell)))]
+global_vertices(mgrid::SimpleMesh3D, cell) = [mgrid.mvertices[v] for v ∈ vertices(cell)]
+
+num_nodes(mgrid::SimpleMesh3D) = length(mgrid.grid.nodes)
+num_faces(mgrid::SimpleMesh3D) = length(mgrid.mfaces)
+num_edges(mgrid::SimpleMesh3D) = length(mgrid.medges)
+num_vertices(mgrid::SimpleMesh3D) = length(mgrid.mvertices)
+
+elementtypes(::SimpleMesh3D{Tetrahedron}) = @SVector [Tetrahedron]
+elementtypes(::SimpleMesh3D{Hexahedron}) = @SVector [Hexahedron]
+
+function to_mesh(grid::Grid{3,C,T}) where {C, T}
+mfaces = OrderedDict{NTuple{3,Int}, Int}()
+medges = OrderedDict{NTuple{2,Int}, Int}()
+mvertices = OrderedDict{Int, Int}()
+next_face_idx = 1
+next_edge_idx = 1
+next_vertex_idx = 1
+number_of_cells_by_type = Dict{DataType, Int}()
+for cell ∈ getcells(grid)
+    cell_type = typeof(cell)
+    if haskey(number_of_cells_by_type, cell_type)
+        number_of_cells_by_type[cell_type] += 1
+    else
+        number_of_cells_by_type[cell_type] = 1
+    end
+
+    for v ∈ vertices(cell)
+        if !haskey(mvertices, v)
+            mvertices[v] = next_vertex_idx
+            next_vertex_idx += 1
+        end
+    end
+    for e ∈ first.(sortedge.(edges(cell)))
+        if !haskey(medges, e)
+            medges[e] = next_edge_idx
+            next_edge_idx += 1
+        end
+    end
+    for f ∈ first.(sortface.(faces(cell)))
+        if !haskey(mfaces, f)
+            mfaces[f] = next_face_idx
+            next_face_idx += 1
+        end
+    end
+end
+return SimpleMesh3D(grid, mfaces, medges, mvertices, number_of_cells_by_type)
+end
+
+
+"""
+SimpleMesh2D{C <: AbstractCell, T <: Real} <: AbstractGrid{2}
+
+A grid which also has information abouts its vertices, faces and edges.
+"""
+struct SimpleMesh2D{C <: AbstractCell, T <: Real} <: AbstractGrid{2}
+grid::Grid{2, C, T}
+mfaces::OrderedDict{NTuple{2,Int}, Int} # Maps "sortface"-representation to id
+mvertices::OrderedDict{Int, Int} # Maps node to id
+number_of_cells_by_type::Dict{DataType, Int}
+end
+
+function to_mesh(grid::Grid{2,C,T}) where {C, T}
+mfaces = OrderedDict{NTuple{2,Int}, Int}()
+mvertices = OrderedDict{Int, Int}()
+next_face_idx = 1
+next_vertex_idx = 1
+number_of_cells_by_type = Dict{DataType, Int}()
+for cell ∈ getcells(grid)
+    cell_type = typeof(cell)
+    if haskey(number_of_cells_by_type, cell_type)
+        number_of_cells_by_type[cell_type] += 1
+    else
+        number_of_cells_by_type[cell_type] = 1
+    end
+
+    for v ∈ vertices(cell)
+        if !haskey(mvertices, v)
+            mvertices[v] = next_vertex_idx
+            next_vertex_idx += 1
+        end
+    end
+    for f ∈ first.(sortface.(faces(cell)))
+        if !haskey(mfaces, f)
+            mfaces[f] = next_face_idx
+            next_face_idx += 1
+        end
+    end
+end
+return SimpleMesh2D(grid, mfaces, mvertices, number_of_cells_by_type)
+end
+
+
+global_faces(mgrid::SimpleMesh2D, cell) = [mgrid.mfaces[sface] for sface ∈ first.(sortface.(faces(cell)))]
+global_vertices(mgrid::SimpleMesh2D, cell) = [mgrid.mvertices[v] for v ∈ vertices(cell)]
+
+num_nodes(mgrid::SimpleMesh2D) = length(mgrid.grid.nodes)
+num_faces(mgrid::SimpleMesh2D) = length(mgrid.mfaces)
+num_vertices(mgrid::SimpleMesh2D) = length(mgrid.mvertices)
+
+elementtypes(::SimpleMesh2D{Triangle}) = @SVector [Triangle]
+elementtypes(::SimpleMesh2D{Quadrilateral}) = @SVector [Quadrilateral]
+
+# Ferrite compat layer for the mesh
+@inline Ferrite.getncells(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getncells(mesh.grid)
+@inline Ferrite.getcells(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getcells(mesh.grid)
+@inline Ferrite.getcells(mesh::Union{SimpleMesh2D,SimpleMesh3D}, v::Union{Int, Vector{Int}}) = Ferrite.getcells(mesh.grid, v)
+@inline Ferrite.getcells(mesh::Union{SimpleMesh2D,SimpleMesh3D}, setname::String) = Ferrite.getcells(mesh.grid, setname)
+@inline Ferrite.getcelltype(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getcelltype(mesh.grid)
+@inline Ferrite.getcelltype(mesh::Union{SimpleMesh2D,SimpleMesh3D}, i::Int) = Ferrite.getcelltype(mesh.grid)
+
+@inline Ferrite.getnnodes(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getnnodes(mesh.grid)
+@inline Ferrite.nnodes_per_cell(mesh::Union{SimpleMesh2D,SimpleMesh3D}, i::Int=1) = Ferrite.nnodes_per_cell(mesh.grid, i)
+@inline Ferrite.getnodes(mesh::Union{SimpleMesh2D,SimpleMesh3D}) = Ferrite.getnodes(mesh.grid)
+@inline Ferrite.getnodes(mesh::Union{SimpleMesh2D,SimpleMesh3D}, v::Union{Int, Vector{Int}}) = Ferrite.getnodes(mesh.grid, v)
+@inline Ferrite.getnodes(mesh::Union{SimpleMesh2D,SimpleMesh3D}, setname::String) = Ferrite.getnodes(mesh.grid, setname)
+
+@inline Ferrite.vtk_grid(filename::AbstractString, mesh::Union{SimpleMesh2D,SimpleMesh3D}; kwargs...) = Ferrite.vtk_grid(filename, mesh.grid, kwargs...)
+
+@inline Ferrite.get_coordinate_type(::SimpleMesh2D{C,T}) where {C,T} = Vec{2,T} 
+@inline Ferrite.get_coordinate_type(::SimpleMesh3D{C,T}) where {C,T} = Vec{3,T} 
diff --git a/src/modeling/common.jl b/src/modeling/common.jl
new file mode 100644
index 000000000..fe270475e
--- /dev/null
+++ b/src/modeling/common.jl
@@ -0,0 +1,11 @@
+# Common modeling primitives are found here
+
+# TODO actually use this :)
+abstract type SteadyStateInternalVariable end
+
+include("core/coefficients.jl")
+
+include("core/boundary_conditions.jl")
+
+include("core/mass.jl")
+include("core/diffusion.jl")
diff --git a/src/modeling/boundary_conditions.jl b/src/modeling/core/boundary_conditions.jl
similarity index 100%
rename from src/modeling/boundary_conditions.jl
rename to src/modeling/core/boundary_conditions.jl
diff --git a/src/modeling/coefficients.jl b/src/modeling/core/coefficients.jl
similarity index 100%
rename from src/modeling/coefficients.jl
rename to src/modeling/core/coefficients.jl
diff --git a/src/modeling/core/diffusion.jl b/src/modeling/core/diffusion.jl
new file mode 100644
index 000000000..4bd0961d5
--- /dev/null
+++ b/src/modeling/core/diffusion.jl
@@ -0,0 +1,38 @@
+# @doc raw"""
+#     AssembledDiffusionOperator{MT, DT, CV}
+
+# Assembles the matrix associated to the bilinearform ``a(u,v) = -\int \nabla v(x) \cdot D(x) \nabla u(x) dx`` for a given diffusion tensor ``D(x)`` and ``u,v`` from the same function space.
+# """
+"""
+    Represents the integrand of the bilinear form <ϕ,ψ> = -∫ D∇ϕ ⋅ ∇ψ dΩ .
+"""
+struct BilinearDiffusionIntegrator{CoefficientType}
+    D::CoefficientType
+    # coordinate_system
+end
+
+struct BilinearDiffusionElementCache{IT <: BilinearDiffusionIntegrator, CV}
+    integrator::IT
+    cellvalues::CV
+end
+
+function assemble_element!(Kₑ, cell, element_cache::CACHE, time) where {CACHE <: BilinearDiffusionElementCache}
+    @unpack cellvalues = element_cache
+    n_basefuncs = getnbasefunctions(cellvalues)
+
+    reinit!(cellvalues, cell)
+
+    for q_point in 1:getnquadpoints(cellvalues)
+        ξ = cellvalues.qr.points[q_point]
+        qp = QuadraturePoint(q_point, ξ)
+        D_loc = evaluate_coefficient(element_cache.integrator.D, cell, qp, time)
+        dΩ = getdetJdV(cellvalues, q_point)
+        for i in 1:n_basefuncs
+            ∇Nᵢ = shape_gradient(cellvalues, q_point, i)
+            for j in 1:n_basefuncs
+                ∇Nⱼ = shape_gradient(cellvalues, q_point, j)
+                Kₑ[i,j] -= ((D_loc ⋅ ∇Nᵢ) ⋅ ∇Nⱼ) * dΩ
+            end
+        end
+    end
+end
diff --git a/src/modeling/core/mass.jl b/src/modeling/core/mass.jl
new file mode 100644
index 000000000..9b57d6a83
--- /dev/null
+++ b/src/modeling/core/mass.jl
@@ -0,0 +1,37 @@
+# @doc raw"""
+#     AssembledMassOperator{MT, CV}
+
+# Assembles the matrix associated to the bilinearform ``a(u,v) = -\int v(x) u(x) dx`` for ``u,v`` from the same function space.
+# """
+"""
+    Represents the integrand of the bilinear form <ϕ,ψ> = ∫ ρϕ ⋅ ψ dΩ .
+"""
+struct BilinearMassIntegrator{CoefficientType}
+    ρ::CoefficientType
+    # coordinate_system
+end
+
+struct BilinearMassElementCache{IT <: BilinearMassIntegrator, T, CV}
+    integrator::IT
+    ρq::Vector{T}
+    cellvalues::CV
+end
+
+function assemble_element!(Mₑ, cell, element_cache::CACHE, time) where {CACHE <: BilinearMassElementCache}
+    @unpack cellvalues = element_cache
+    reinit!(element_cache.cellvalues, cell)
+    n_basefuncs = getnbasefunctions(cellvalues)
+    for q_point in 1:getnquadpoints(cellvalues)
+        ξ = cellvalues.qr.points[q_point]
+        qp = QuadraturePoint(q_point, ξ)
+        ρ = evaluate_coefficient(element_cache.integrator.ρ, cell, qp, time)
+        dΩ = getdetJdV(cellvalues, q_point)
+        for i in 1:n_basefuncs
+            Nᵢ = shape_value(cellvalues, q_point, i)
+            for j in 1:n_basefuncs
+                Nⱼ = shape_value(cellvalues, q_point, j)
+                Mₑ[i,j] += ρ * Nᵢ * Nⱼ * dΩ 
+            end
+        end
+    end
+end
diff --git a/src/modeling/electrophysiology.jl b/src/modeling/electrophysiology.jl
index 8a36f565e..d3eee3a66 100644
--- a/src/modeling/electrophysiology.jl
+++ b/src/modeling/electrophysiology.jl
@@ -143,6 +143,7 @@ The original model formulation (TODO citation) with the structure
     ∂ₜ𝐬  = g(φₘ,𝐬,x)
  φᵢ - φₑ = φₘ
 
+TODO implement
 """
 struct ParabolicParabolicBidomainModel <: AbstractEPModel
     χ
@@ -163,6 +164,8 @@ Transformed bidomain model with the structure
 
 This formulation is a transformation of the parabolic-parabolic
 form (c.f. TODO citation) and has been derived by (TODO citation) first.
+
+TODO implement
 """
 struct ParabolicEllipticBidomainModel <: AbstractEPModel
     χ
@@ -173,7 +176,6 @@ struct ParabolicEllipticBidomainModel <: AbstractEPModel
     ion::AbstractIonicModel
 end
 
-
 """
 Simplification of the bidomain model with the structure
 
@@ -192,84 +194,8 @@ struct MonodomainModel{F1,F2,F3,STIM<:TransmembraneStimulationProtocol,ION<:Abst
     ion::ION
 end
 
-# @doc raw"""
-#     AssembledDiffusionOperator{MT, DT, CV}
-
-# Assembles the matrix associated to the bilinearform ``a(u,v) = -\int \nabla v(x) \cdot D(x) \nabla u(x) dx`` for a given diffusion tensor ``D(x)`` and ``u,v`` from the same function space.
-# """
-"""
-    Represents the integrand of the bilinear form <ϕ,ψ> = -∫ D∇ϕ ⋅ ∇ψ dΩ .
-"""
-struct BilinearDiffusionIntegrator{CoefficientType}
-    D::CoefficientType
-    # coordinate_system
-end
-
-struct BilinearDiffusionElementCache{IT <: BilinearDiffusionIntegrator, CV}
-    integrator::IT
-    cellvalues::CV
-end
-
-function assemble_element!(Kₑ, cell, element_cache::CACHE, time) where {CACHE <: BilinearDiffusionElementCache}
-    @unpack cellvalues = element_cache
-    n_basefuncs = getnbasefunctions(cellvalues)
-
-    reinit!(cellvalues, cell)
-
-    for q_point in 1:getnquadpoints(cellvalues)
-        ξ = cellvalues.qr.points[q_point]
-        qp = QuadraturePoint(q_point, ξ)
-        D_loc = evaluate_coefficient(element_cache.integrator.D, cell, qp, time)
-        dΩ = getdetJdV(cellvalues, q_point)
-        for i in 1:n_basefuncs
-            ∇Nᵢ = shape_gradient(cellvalues, q_point, i)
-            for j in 1:n_basefuncs
-                ∇Nⱼ = shape_gradient(cellvalues, q_point, j)
-                Kₑ[i,j] -= ((D_loc ⋅ ∇Nᵢ) ⋅ ∇Nⱼ) * dΩ
-            end
-        end
-    end
-end
-
-# @doc raw"""
-#     AssembledMassOperator{MT, CV}
-
-# Assembles the matrix associated to the bilinearform ``a(u,v) = -\int v(x) u(x) dx`` for ``u,v`` from the same function space.
-# """
-"""
-    Represents the integrand of the bilinear form <ϕ,ψ> = ∫ ρϕ ⋅ ψ dΩ .
-"""
-struct BilinearMassIntegrator{CoefficientType}
-    ρ::CoefficientType
-    # coordinate_system
-end
-
-struct BilinearMassElementCache{IT <: BilinearMassIntegrator, T, CV}
-    integrator::IT
-    ρq::Vector{T}
-    cellvalues::CV
-end
-
-function assemble_element!(Mₑ, cell, element_cache::CACHE, time) where {CACHE <: BilinearMassElementCache}
-    @unpack cellvalues = element_cache
-    reinit!(element_cache.cellvalues, cell)
-    n_basefuncs = getnbasefunctions(cellvalues)
-    for q_point in 1:getnquadpoints(cellvalues)
-        ξ = cellvalues.qr.points[q_point]
-        qp = QuadraturePoint(q_point, ξ)
-        ρ = evaluate_coefficient(element_cache.integrator.ρ, cell, qp, time)
-        dΩ = getdetJdV(cellvalues, q_point)
-        for i in 1:n_basefuncs
-            Nᵢ = shape_value(cellvalues, q_point, i)
-            for j in 1:n_basefuncs
-                Nⱼ = shape_value(cellvalues, q_point, j)
-                Mₑ[i,j] += ρ * Nᵢ * Nⱼ * dΩ 
-            end
-        end
-    end
-end
-
 ######################################################
+# TODO where to put these?
 Base.@kwdef struct ParametrizedFHNModel{T} <: AbstractIonicModel
     a::T = T(0.1)
     b::T = T(0.5)
diff --git a/src/modeling/fluid_mechanics.jl b/src/modeling/fluid_mechanics.jl
new file mode 100644
index 000000000..04f970ec0
--- /dev/null
+++ b/src/modeling/fluid_mechanics.jl
@@ -0,0 +1 @@
+include("fluid/lumped.jl")
diff --git a/src/modeling/multiphysics.jl b/src/modeling/multiphysics.jl
new file mode 100644
index 000000000..095af7a8c
--- /dev/null
+++ b/src/modeling/multiphysics.jl
@@ -0,0 +1,2 @@
+include("coupler/interface.jl")
+include("coupler/fsi.jl")
diff --git a/src/modeling/solid/elements.jl b/src/modeling/solid/elements.jl
new file mode 100644
index 000000000..51da8e5e1
--- /dev/null
+++ b/src/modeling/solid/elements.jl
@@ -0,0 +1,61 @@
+
+"""
+TODO document me!
+"""
+# Is this really what we want? Or should the face integrals be directly merged into the mechanical model?
+# Does this abstraction layer have another benefit?
+struct StructuralModel{MM, FM}
+    mechanical_model::MM
+    face_models::FM
+end
+
+
+"""
+    StructuralElementCache
+
+TODO document me!
+"""
+# TODO rename contraction_model_cache
+struct StructuralElementCache{M, CMCache, CV}
+    constitutive_model::M
+    contraction_model_cache::CMCache
+    cv::CV
+end
+
+# TODO how to control dispatch on required input for the material routin?
+# TODO finer granularity on the dispatch here. depending on the evolution law of the internal variable this routine looks slightly different.
+function assemble_element!(Kₑ::Matrix, residualₑ, uₑ, geometry_cache, element_cache::StructuralElementCache, time)
+    @unpack constitutive_model, contraction_model_cache, cv = element_cache
+    ndofs = getnbasefunctions(cv)
+
+    reinit!(cv, geometry_cache)
+
+    @inbounds for qpᵢ in 1:getnquadpoints(cv)
+        ξ = cv.qr.points[qpᵢ]
+        qp = QuadraturePoint(qpᵢ, ξ)
+        dΩ = getdetJdV(cv, qpᵢ)
+
+        # Compute deformation gradient F
+        ∇u = function_gradient(cv, qpᵢ, uₑ)
+        F = one(∇u) + ∇u
+
+        # Compute stress and tangent
+        contraction_state = state(contraction_model_cache, geometry_cache, qp, time)
+        P, ∂P∂F = material_routine(constitutive_model, F, contraction_state, geometry_cache, qp, time)
+
+        # Loop over test functions
+        for i in 1:ndofs
+            ∇δui = shape_gradient(cv, qpᵢ, i)
+
+            # Add contribution to the residual from this test function
+            residualₑ[i] += ∇δui ⊡ P * dΩ
+
+            ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation
+            for j in 1:ndofs
+                ∇δuj = shape_gradient(cv, qpᵢ, j)
+                # Add contribution to the tangent
+                Kₑ[i, j] += ( ∇δui∂P∂F ⊡ ∇δuj ) * dΩ
+            end
+        end
+    end
+end
diff --git a/src/modeling/solid/drivers.jl b/src/modeling/solid/materials.jl
similarity index 52%
rename from src/modeling/solid/drivers.jl
rename to src/modeling/solid/materials.jl
index efca324ec..0bfd746c6 100644
--- a/src/modeling/solid/drivers.jl
+++ b/src/modeling/solid/materials.jl
@@ -1,17 +1,14 @@
-
-struct StructuralModel{MM, FM}
-    mechanical_model::MM
-    face_models::FM
-end
+# TODO (FILE) I think we should change the design here. Instea of dispatching on Ψ we should make the material callable or equip it with a function.
 
 abstract type QuasiStaticModel end
 
-#TODO constrain to orthotropic material models, e.g. via traits, or rewrite all 3 "constitutive_driver"s below
-function constitutive_driver(constitutive_model, F, internal_state, geometry_cache::Ferrite.CellCache, qp::QuadraturePoint, time)
+#TODO constrain to orthotropic material models, e.g. via traits, or rewrite all 3 "material_routine"s below
+function material_routine(constitutive_model, F, internal_state, geometry_cache::Ferrite.CellCache, qp::QuadraturePoint, time)
     f₀, s₀, n₀ = evaluate_coefficient(constitutive_model.microstructure_model, geometry_cache, qp, time)
-    return constitutive_driver(F, f₀, s₀, n₀, internal_state, constitutive_model)
+    return material_routine(F, f₀, s₀, n₀, internal_state, constitutive_model)
 end
 
+
 """
 @TODO citation
 """
@@ -25,7 +22,7 @@ end
 
 """
 """
-function constitutive_driver(F::Tensor{2,dim}, f₀::Vec{dim}, s₀::Vec{dim}, n₀::Vec{dim}, internal_state, model::GeneralizedHillModel) where {dim}
+function material_routine(F::Tensor{2,dim}, f₀::Vec{dim}, s₀::Vec{dim}, n₀::Vec{dim}, internal_state, model::GeneralizedHillModel) where {dim}
     # TODO what is a good abstraction here?
     Fᵃ = compute_Fᵃ(internal_state, f₀, s₀, n₀, model.contraction_model, model.active_deformation_gradient_model)
 
@@ -38,6 +35,7 @@ function constitutive_driver(F::Tensor{2,dim}, f₀::Vec{dim}, s₀::Vec{dim}, n
     return ∂Ψ∂F, ∂²Ψ∂F²
 end
 
+
 """
 @TODO citation
 """
@@ -51,7 +49,7 @@ end
 
 """
 """
-function constitutive_driver(F::Tensor{2,dim}, f₀::Vec{dim}, s₀::Vec{dim}, n₀::Vec{dim}, cell_state, model::ExtendedHillModel) where {dim}
+function material_routine(F::Tensor{2,dim}, f₀::Vec{dim}, s₀::Vec{dim}, n₀::Vec{dim}, cell_state, model::ExtendedHillModel) where {dim}
     # TODO what is a good abstraction here?
     Fᵃ = compute_Fᵃ(cell_state, f₀, s₀, n₀, model.contraction_model, model.active_deformation_gradient_model)
     N = 𝓝(cell_state, model.contraction_model)
@@ -65,6 +63,7 @@ function constitutive_driver(F::Tensor{2,dim}, f₀::Vec{dim}, s₀::Vec{dim}, n
     return ∂Ψ∂F, ∂²Ψ∂F²
 end
 
+
 """
 """
 struct ActiveStressModel{Mat, ASMod, CMod, MS} <: QuasiStaticModel
@@ -76,7 +75,7 @@ end
 
 """
 """
-function constitutive_driver(F::Tensor{2,dim}, f₀::Vec{dim}, s₀::Vec{dim}, n₀::Vec{dim}, cell_state, model::ActiveStressModel) where {dim}
+function material_routine(F::Tensor{2,dim}, f₀::Vec{dim}, s₀::Vec{dim}, n₀::Vec{dim}, cell_state, model::ActiveStressModel) where {dim}
     ∂²Ψ∂F², ∂Ψ∂F = Tensors.hessian(
         F_ad ->
               Ψ(F_ad,     f₀, s₀, n₀, model.material_model),
@@ -97,47 +96,3 @@ struct ElastodynamicsModel{RHSModel <: QuasiStaticModel, CoefficienType}
     rhs::RHSModel
     ρ::CoefficienType
 end
-
-"""
-"""
-struct StructuralElementCache{M, CMCache, CV}
-    constitutive_model::M
-    contraction_model_cache::CMCache
-    cv::CV
-end
-
-function assemble_element!(Kₑ::Matrix, residualₑ, uₑ, geometry_cache, element_cache::StructuralElementCache, time)
-    @unpack constitutive_model, contraction_model_cache, cv = element_cache
-    ndofs = getnbasefunctions(cv)
-
-    reinit!(cv, geometry_cache)
-
-    @inbounds for qpᵢ in 1:getnquadpoints(cv)
-        ξ = cv.qr.points[qpᵢ]
-        qp = QuadraturePoint(qpᵢ, ξ)
-        dΩ = getdetJdV(cv, qpᵢ)
-
-        # Compute deformation gradient F
-        ∇u = function_gradient(cv, qpᵢ, uₑ)
-        F = one(∇u) + ∇u
-
-        # Compute stress and tangent
-        contraction_state = state(contraction_model_cache, geometry_cache, qp, time)
-        P, ∂P∂F = constitutive_driver(constitutive_model, F, contraction_state, geometry_cache, qp, time)
-
-        # Loop over test functions
-        for i in 1:ndofs
-            ∇δui = shape_gradient(cv, qpᵢ, i)
-
-            # Add contribution to the residual from this test function
-            residualₑ[i] += ∇δui ⊡ P * dΩ
-
-            ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation
-            for j in 1:ndofs
-                ∇δuj = shape_gradient(cv, qpᵢ, j)
-                # Add contribution to the tangent
-                Kₑ[i, j] += ( ∇δui∂P∂F ⊡ ∇δuj ) * dΩ
-            end
-        end
-    end
-end
diff --git a/src/modeling/solid_mechanics.jl b/src/modeling/solid_mechanics.jl
new file mode 100644
index 000000000..01877f6cd
--- /dev/null
+++ b/src/modeling/solid_mechanics.jl
@@ -0,0 +1,5 @@
+include("solid/energies.jl")
+include("solid/contraction.jl")
+include("solid/active.jl")
+include("solid/materials.jl")
+include("solid/elements.jl")
diff --git a/src/solver/backward_euler.jl b/src/solver/euler.jl
similarity index 100%
rename from src/solver/backward_euler.jl
rename to src/solver/euler.jl
diff --git a/src/utils.jl b/src/utils.jl
index 1ce021bd6..aee88b95c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,3 +1,6 @@
+
+include("collections.jl")
+
 function calculate_element_volume(cell, cellvalues_u, uₑ)
     reinit!(cellvalues_u, cell)
     evol::Float64=0.0;
diff --git a/test/test_type_stability.jl b/test/test_type_stability.jl
index b971540e9..f931fa35e 100644
--- a/test/test_type_stability.jl
+++ b/test/test_type_stability.jl
@@ -33,7 +33,7 @@
                     PelceSunLangeveld1995Model(;calcium_field=ConstantCoefficient(1.0)),
                     fsn,
                 )
-                @test_call Thunderbolt.constitutive_driver(F, f₀, s₀, n₀, Caᵢ, model)
+                @test_call Thunderbolt.material_routine(F, f₀, s₀, n₀, Caᵢ, model)
             end
         end
 
@@ -57,7 +57,7 @@
                             contraction_model,
                             fsn,
                         )
-                        @test_call broken=true Thunderbolt.constitutive_driver(F, f₀, s₀, n₀, Caᵢ, model)
+                        @test_call broken=true Thunderbolt.material_routine(F, f₀, s₀, n₀, Caᵢ, model)
                     end
                 end
             end