Skip to main content

A convenient scheme for coupling a finite element curvilinear mesh to a finite element voxel mesh: application to the heart



In some cases, it may be necessary to combine distinct finite element meshes into a single system. The present work describes a scheme for coupling a finite element mesh, which may have curvilinear elements, to a voxel based finite element mesh.


The method is described with reference to a sample problem that involves combining a heart, which is defined by a curvilinear mesh, with a voxel based torso mesh. The method involves the creation of a temporary (scaffolding) mesh that couples the outer surface of the heart mesh to a voxel based torso mesh. The inner surface of the scaffolding mesh is the outer heart surface, and the outer surface of the scaffolding mesh is defined by the nodes in the torso mesh that are nearest (but outside of) the heart. The finite element stiffness matrix for the scaffolding mesh is then computed. This stiffness matrix includes extraneous nodes that are then removed, leaving a coupling matrix that couples the original outer heart surface nodes to adjacent nodes in the torso voxel mesh. Finally, a complete system matrix is assembled from the pre-existing heart stiffness matrix, the heart/torso coupling matrix, and the torso stiffness matrix.


Realistic body surface electrocardiograms were generated. In a test involving a dipole embedded in a spherical shell, relative error of the scheme rapidly converged to slightly over 4%, although convergence thereafter was relatively slow.


The described method produces reasonably accurate results and may be best suited for problems where computational speed and convenience have a higher priority than very high levels of accuracy.


In some cases, it may be necessary to combine distinct finite element meshes into a single system. For example, the task that motivated this work involved integrating a curvilinear mesh[1] of the Auckland canine heart[2] with a voxel mesh of a human shaped torso. The heart mesh described in [1] was created to perform simulations on a stand alone heart. Fischer et al. [3] describe an elegant technique for coupling a stand alone heart model to a torso surface via the boundary element method, but this scheme is not applicable where it is desired to know the intra-thoracic potential distribution, which requires a full volume mesh.

Many types of volume meshing strategies exist[4]. Meshes can generally be divided into two categories, structured meshes and unstructured meshes. Structured meshes, which include voxel meshes, are characterized by interior nodes that are connected to the same number of nodes. In contrast, the internal nodes in an unstructured mesh may be connected to different numbers of nodes. Structured meshes are generally easy to work with because of their regularity. However, it may be difficult to conform a structured mesh to an irregular boundary, a task which is easier with an unstructured mesh.

One approach for coupling a voxel mesh to a mesh with an irregular boundary (e.g. the outer heart surface) involves individual treatment of the voxel elements that intersect the boundary. This scheme, known as the Cartesian cut cell method[5], has the drawback of being relatively complex and therefore likely to be difficult to implement. Another coupling technique involves extruding prism elements from a triangular surface mesh (e.g. the outer heart surface) in the surface normal direction, and then coupling the prism elements to Cartesian elements[6]. Again, the disadvantage of this approach is that it is relatively complicated.

The present work describes a technique that is a variant of the hybrid Cartesian/prism scheme, with a novel method for coupling the prism elements to the voxel elements. The scheme is easy to code, computes quickly, requires little or no user supervision, and has all the advantages associated with voxel meshes. The main drawbacks of the scheme are relatively uneven and slow convergence after a certain degree of accuracy has been reached. One additional limitation is that the scheme is most suitable for geometries which are reasonably well characterized by a spherical harmonic expansion.



The present scheme will be described with reference to a particular problem that involves coupling a previously meshed heart[1] to a torso, as shown in the left hand panel in Figure 1. The goal is to solve for the extracellular potential (V e) throughout the heart and the potential outside of the heart (V o) in the elliptical equation[7]

Figure 1
figure 1

Sample problem. The heart/torso system is shown in the left panel. The right panel shows corresponding simulated electrograms at the sites designated in the left panel.

•(D i + D e )V e = -•(D i )V t ; within heart muscle

•(σ o )V o = 0; in the torso and blood pools     (1)

where D i and D e are the intracellular and extracellular conductivity tensors within heart muscle, σo represents the conductivity in the torso and blood pools respectively, and V t are the transmembrane potentials throughout heart muscle. At the interface between the heart and torso, the potential is continuous (Ve = Vo), the normal component of the extracellular current is continuous with the normal component of the current through the torso/blood pools (nT D eVe = nT σo Vo, where n is the unit normal vector), and there is no flow of intracellular current into the torso/blood pools (nT D iVi = 0).

Figure 2 is a flow chart of the coupling scheme. The first step of the method involves embedding the heart in a preliminary voxel mesh that is slightly larger than the heart, as indicated in Figure 3. The next step involves locating preliminary voxel mesh nodes ("boundary nodes") that are outside of the heart but close to the heart, as indicated by the red circles in Figure 3. A surface ("boundary surface") is fitted to these boundary nodes. The boundary surface is essentially a projection of the outer heart surface, as shown in Figure 3. A temporary mesh ("scaffolding mesh") is generated between the outer heart surface and the boundary surface and the stiffness matrix for this mesh is computed. The scaffolding mesh is indicated by the purple area in Figure 3. The scaffolding mesh includes extraneous temporary nodes that are then removed from the stiffness matrix, leaving a coupling matrix that couples the pre-existing outer heart surface nodes to the boundary surface voxel nodes. A voxel mesh is generated for the entire torso by adding voxels to the preliminary mesh. The resulting stiffness matrix for the torso nodes is then computed. Finally, a complete system matrix is assembled from the pre-existing heart stiffness matrix, the heart/torso coupling matrix and the torso stiffness matrix.

Figure 2
figure 2

Flow chart of the coupling scheme.

Figure 3
figure 3

General mesh structure in 2 dimensions. The heart is embedded in a preliminary "scaffolding" voxel mesh, which comprises all of the voxels within the blue rectangle. The preliminary mesh is in turn embedded and aligned with a voxel torso mesh, whose outer nodes are shown in green. Within the preliminary mesh, the voxels that have nodes both inside and outside of the heart are located, and the nodes in these voxels that are outside of the heart are defined as boundary nodes (Xb). A surface (the "boundary surface") is fitted to these nodes. A coupling mesh (in purple) between the outer heart surface and the boundary surface is then generated.

Preliminary mesh; boundary extraction

Given the heart surface, a preliminary voxel mesh is generated that is slightly larger than the heart. In the two dimensional example illustrated in Figure 3, the preliminary voxel mesh encloses all of the squares within the blue rectangle. The nodes of this mesh are designated as Xp = (xp, yp, zp). It will be assumed that the voxels are isotropic (side length = s) and chosen to be defined with respect to some convenient origin. The preliminary voxel mesh should extend for a distance of at least 1s in the x, y and z directions beyond any point on the outer heart surface. The voxel size s is preferably the desired voxel size for the entire torso mesh.

The next step involves locating boundary nodes Xb Xp that are outside of the heart but closest to it, as illustrated in Figures 3 and 4. The boundary nodes are located by finding boundary voxels Vb, which are defined as those voxels that have at least one node within the heart and at least one node outside of the heart. The boundary nodes Xb Xp are defined as all of the nodes within these voxels Vb that are outside of the heart.

Figure 4
figure 4

Finding boundary nodes and creating the boundary surface. A convenient origin within the heart serves as the basis for fitting a spherical harmonic expansion to the outer heart surface, which generates a surface characterized by the radius Rh(θ,Φ). Voxels that have nodes both within and without of this surface are boundary voxels. Within these voxels, the nodes outside of the heart are boundary nodes (red). A surface ("boundary surface") is fitted to these nodes by performing another spherical harmonic expansion about the heart origin. This surface is characterized by the radius Rb(θ,Φ).

There are a number of methods for locating the boundary voxels. One convenient method, which was used in the present example, is applicable to any surface (e.g. the heart surface) that may be characterized by spherical harmonics[1]. In this case, the spherical harmonics provide a continuous representation of a surface, and therefore may be used to determine whether any given point is inside or outside of this surface.

More particularly, the outer heart surface Sh is approximated by Rh = F(θ,Φ), where Rh is the radius of the approximate outer heart surface, F is a spherical harmonic polynomial as will be further described below, and θ and Φ are azimuthal and polar angles, respectively. Figure 4 shows a 2-D representation of Rh. Rh, θ, and Φ are defined with respect to a convenient origin which is at the point (xh, yh, zh). This origin is preferably chosen such that it results in a good spherical harmonic expansion. The grid points Xp are translated so their origin is (xh, yh, zh), which results in a set of translated node points X'p = (xp- xh, yp - yh, zp - zh). The azimuthal and polar angles θp and Φp (with respect to the (xh, yh, zh) origin) are computed for X'p, and then the value F(θp, Φp) is computed, which is the radius of the outer heart surface at the angles θp and Φp. If this radius is greater than the radius ( x p x h ) 2 + ( y p y h ) 2 + ( z p z h ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaadaGcaaqaaiabcIcaOiabdIha4naaBaaaleaacqWGWbaCaeqaaOGaeyOeI0IaemiEaG3aaSbaaSqaaiabdIgaObqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaakiabgUcaRiabcIcaOiabdMha5naaBaaaleaacqWGWbaCaeqaaOGaeyOeI0IaemyEaK3aaSbaaSqaaiabdIgaObqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaakiabgUcaRiabcIcaOiabdQha6naaBaaaleaacqWGWbaCaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdIgaObqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeqaaaaa@4C3A@ of a translated node point, then that point is within the heart. The set of node points within the heart will be denoted by the set Xh Xp.

The above procedure locates the nodes Xh Xp within the heart. The boundary voxels Vb are then defined as those voxels that have at least one node within the heart (Xh) and at least one node outside of the heart. The boundary points Xb are those nodes that are within the boundary voxels but outside of the heart.

Surface formation

A surface (the boundary surface Sb) is generated by fitting a spherical harmonic expansion to the points Xb about the origin (xh, yh, zh). The radii (Rb) and azimuthal and polar angles θb and Φb are computed for the boundary points, as indicated in Figure 4, which shows only the azimuthal angle. The spherical harmonic expansion of the radii Rb may be expressed as

R b = i = 0 N { a i P i ( cos ϕ b ) + m = 1 i ( a i m cos ( m θ b ) + b i m sin ( m θ b ) ) P i m ( cos ϕ b ) } R b = P C ( 2 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaafaqaaeGabaaabaGaemOuai1aaSbaaSqaaiabdkgaIbqabaGccqGH9aqpdaaeWbqaamaacmqabaGaemyyae2aaSbaaSqaaiabdMgaPbqabaGccqWGqbaudaWgaaWcbaGaemyAaKgabeaakiabcIcaOiGbcogaJjabc+gaVjabcohaZHGaciab=v9aQnaaBaaaleaacqWGIbGyaeqaaOGaeiykaKIaey4kaSYaaabCaeaacqGGOaakcqWGHbqydaWgaaWcbaGaemyAaKMaemyBa0gabeaakiGbcogaJjabc+gaVjabcohaZjabcIcaOiabd2gaTjab=H7aXnaaBaaaleaacqWGIbGyaeqaaOGaeiykaKIaey4kaSIaemOyai2aaSbaaSqaaiabdMgaPjabd2gaTbqabaGccyGGZbWCcqGGPbqAcqGGUbGBcqGGOaakcqWGTbqBcqWF4oqCdaWgaaWcbaGaemOyaigabeaakiabcMcaPiabcMcaPiabdcfaqnaaDaaaleaacqWGPbqAaeaacqWGTbqBaaGccqGGOaakcyGGJbWycqGGVbWBcqGGZbWCcqWFvpGAdaWgaaWcbaGaemOyaigabeaakiabcMcaPaWcbaGaemyBa0Maeyypa0JaeGymaedabaGaemyAaKganiabggHiLdaakiaawUhacaGL9baaaSqaaiabdMgaPjabg2da9iabicdaWaqaaiabd6eaobqdcqGHris5aaGcbaacbeGae4Nuai1aaSbaaSqaaiab+jgaIbqabaGccqGH9aqpcqGFqbaucqGHxiIkcqGFdbWqaaGaaCzcaiaaxMaadaqadaqaaiabikdaYaGaayjkaiaawMcaaaaa@87CC@

where N is the number of terms in the spherical harmonic expansion, Pi is the Legendre function of degree i of the first kind, Pi m is the associated Legendre function of degree i and order m, ai, aim, bim are coefficients (to be determined), P and C are matrix and vector representations of the Legendre function values and coefficients, respectively, and R b is a vector representation of the radii Rb. A least squares fit (C = P\R b in Matlab notation) results in a set of coefficients Cbs that define the boundary surface Sb characterized by Rbs = F(θ,Φ), where F is the spherical harmonic expansion of order N with coefficients Cbs, and Rbs is the radius of a point on the surface Sb at a given pair of azimuthal and polar angles (θ, Φ).

The appropriate value of N will depend on the geometry of the problem and may be determined empirically. For the heart shown in Figure 1, N = 8 produced good results.

Scaffolding mesh generation

The next step involves forming a scaffolding mesh that couples the original heart nodes Xh to the boundary nodes Xb. The scaffolding mesh, which is indicated by the purple region in Figure 3, joins the heart surface Sh to the boundary surface Sb. With reference to Figure 5, which again shows a two dimensional example of the present scheme, the scaffolding mesh is shown in light purple and one of the scaffolding mesh elements is indicated by a dotted black quadrilateral. The scaffolding mesh provides an approximation of the electrical coupling between the original heart mesh nodes (Xh, shown as black filled circles) and the boundary nodes (Xb, shown as red filled circles). The scaffolding mesh is only an estimate of this electrical coupling because the mesh does not include the actual Xb nodes, but approximations (dotted red circles) of the locations of the boundary nodes (Xb).

Figure 5
figure 5

Creating the scaffolding mesh. The boundary nodes (Xb) are projected onto both the heart surface and the boundary surface, as indicated by the red diagonal line circles and red dotted circles, respectively. The original heart surface nodes are projected on to the boundary surface (black diagonal lined circle). The projected nodes along with the original heart surface node define an element of the scaffolding mesh (shown in black dotted fill for just one element).

In three dimensions, the scaffolding mesh comprises prismatic elements. Each prismatic element consists of two similar triangles, one on the heart surface Sh and one on the boundary surface Sb, whose corresponding vertices are joined to one another. The triangles, in turn, are formed from combinations of the original heart nodes Xh and boundary nodes Xb, and projections of these nodes.

More particularly, the heart nodes Xh are projected onto the surface Sb by a line that passes through the heart nodes and the heart origin (xh, yh, zh), as illustrated in Figure 4. Similarly, the boundary nodes Xb are projected onto both the surface Sb and the surface Sh by a line that passes through these nodes and the heart origin (xh, yh, zh). The projection of the boundary nodes on to Sb will be referred to as Xb'. From an implementation standpoint, the projections may be accomplished by calculating the (θ,Φ) pair for any point (about the (xh, yh, zh) origin), and then computing the radius R = F(θ,Φ) for either the heart surface Sh or boundary surface Sb from its corresponding spherical harmonic expansion F(θ,Φ).

The result of the mutual projections is two sets of points, one on Sh and the other on Sb, such that each point at a given (θ,Φ) on a surface has a corresponding point on the other surface at the same (θ,Φ). Thus, a triangular tessellation of one surface may be used to generate an identical triangular tessellation of the other surface.

There are many ways to tessellate a surface. The approach that was adopted took advantage of the regularity of the (θ,Φ) coordinates of the heart nodes Xh. These nodes formed a regular grid in the (θ,Φ) coordinates by choice; it is easy to form a such a regular grid when the surface at issue (Sb) is characterized by spherical harmonics. Because of the cyclicity of the azimuthal angle θ, the original (θ,Φ) grid (-π<=θ<=π) points were augmented such that those points (θa, Φa) with θa<=-4π/5 were identified, and new points with (θa+2*π, Φa) added to the grid. This "wrapped around" grid was triangulated, and the triangles characterized by all three points having θ>2*π were removed. This type of tessellation is generally acceptable only in cases when there are points at both poles (Φ = +/-pi) and there is a reasonably smooth gradation in the Φ values of the points. This condition was satisfied in the present case because the original heart points Xh were generated via a spherical harmonic expansion with an even division in the Φ direction.

The left panel in Figure 6 shows the surface of the example problem heart with both the original heart surface points (blue) and projected grid boundary points Xb (red) triangulated in the manner described above. The boundary surface (right hand panel in Figure) has shape very similar to that of the heart surface and the triangulation is identical to that shown for the heart surface.

Figure 6
figure 6

Tessellated surfaces. The heart surface is in the left panel and the boundary surface is in the right panel. The tessellated patterns for the surfaces are identical. In the right hand panel, the boundary surface is represented by a transparent mesh that shows the embedded heart surface (aqua).

Because each triangle on the heart surface has a corresponding triangle on the boundary surface, connecting the corresponding triangle nodes by radial line segments defines a prismatic element. All of the prismatic elements define the scaffolding mesh.

Heart to torso transfer matrix

The finite element stiffness matrix (Ak) for the prismatic element (scaffolding) mesh describes the electrical coupling between the heart surface and the boundary surface. However, this matrix contains extraneous elements, i.e. the projected heart nodes on the boundary surface (Xhp), and the projected boundary nodes on the heart surface (Xbp). These extraneous nodes must be removed, resulting in a transfer matrix between the heart nodes and the boundary nodes.

The concept of extraneous nodes will be further described with reference to Figure 5. The desired final matrix should reflect only electrical coupling between the original heart mesh nodes (Xh) shown as filled black circles and the boundary nodes (Xb) shown as filled red circles. The positions of the boundary nodes (Xb) are approximated by the dotted red circles (Xb'), so in actuality the final desired stiffness matrix will reflect the coupling between the filled black circles and the dotted red circles. However, the scaffolding matrix also contains electrical coupling between the desired nodes (filled black and dotted red circles) and extraneous nodes (circles with black and red diagonal lines), which are the projected heart nodes on the boundary surface (Xhp), and the projected boundary nodes on the heart surface (Xbp). These extraneous nodes must not exist in the desired final system matrix.

First, to assemble the finite element stiffness matrix (Ak) for the prismatic element mesh, the local stiffness matrix for each prism element is constructed. This may be done in a number of ways. For a Galerkin type FEM, a closed form local matrix based on the prismatic geometry may be used. Alternatively, each prism may be decomposed into tetrahedrons, and the stiffness matrix, based on linear basis functions, may be computed for each tetrahedron.

To remove the extraneous nodes in Ak, it was partitioned as follows:

A k = [ A h h k A h e k A h b k A e h k A e e k A e b k A b h k A b e k A b b k ] ( 3 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaahaaWcbeqaaiabdUgaRbaakiabg2da9maadmaabaqbaeqabmWaaaqaaGqadiab+feabnaaDaaaleaacqWGObaAcqWGObaAaeaacqWGRbWAaaaakeaacqWFbbqqdaqhaaWcbaGaemiAaGMaemyzaugabaGaem4AaSgaaaGcbaGae8xqae0aa0baaSqaaiabdIgaOjqbdkgaIzaafaaabaGaem4AaSgaaaGcbaGae8xqae0aa0baaSqaaiabdwgaLjabdIgaObqaaiabdUgaRbaaaOqaaiab=feabnaaDaaaleaacqWGLbqzcqWGLbqzaeaacqWGRbWAaaaakeaacqWFbbqqdaqhaaWcbaGaemyzauMafmOyaiMbauaaaeaacqWGRbWAaaaakeaacqWFbbqqdaqhaaWcbaGafmOyaiMbauaacqWGObaAaeaacqWGRbWAaaaakeaacqWFbbqqdaqhaaWcbaGafmOyaiMbauaacqWGLbqzaeaacqWGRbWAaaaakeaacqWFbbqqdaqhaaWcbaGafmOyaiMbauaacuWGIbGygaqbaaqaaiabdUgaRbaaaaaakiaawUfacaGLDbaacaWLjaGaaCzcamaabmaabaGaeG4mamdacaGLOaGaayzkaaaaaa@65D1@

where the subscript h denotes the heart surface nodes Xh, the subscript e denotes the extraneous nodes (Xhp Xbp) to be removed, and the subscript b' denotes the projected boundary nodes Xb'. The desired coupling matrix (Ac) between Xh and Xb' is found by solving for the extraneous nodes in the set of equations Ak = 0:

A c = [ A h h c A h b c A b h c A b b c ] ( 4 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaahaaWcbeqaaiabdogaJbaakiabg2da9maadmaabaqbaeqabiGaaaqaaiab=feabnaaDaaaleaacqWGObaAcqWGObaAaeaacqWGJbWyaaaakeaacqWFbbqqdaqhaaWcbaGaemiAaGMafmOyaiMbauaaaeaacqWGJbWyaaaakeaacqWFbbqqdaqhaaWcbaGafmOyaiMbauaacqWGObaAaeaacqWGJbWyaaaakeaacqWFbbqqdaqhaaWcbaGafmOyaiMbauaacuWGIbGygaqbaaqaaiabdogaJbaaaaaakiaawUfacaGLDbaacaWLjaGaaCzcamaabmaabaGaeGinaqdacaGLOaGaayzkaaaaaa@4B0A@


A h h c = A h h k A h e k ( A e e k ) 1 A e h k A h b c = A h b k A h e k ( A e e k ) 1 A e b k A b h c = A b h k A b e k ( A e e k ) 1 A e h k A b b c = A b b k A b e k ( A e e k ) 1 A e b k MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaafaqaaeabbaaaaeaaieqacqWFbbqqdaqhaaWcbaGaemiAaGMaemiAaGgabaGaem4yamgaaOGaeyypa0Jae8xqae0aa0baaSqaaiabdIgaOjabdIgaObqaaiabdUgaRbaakiabgkHiTiab=feabnaaDaaaleaacqWGObaAcqWGLbqzaeaacqWGRbWAaaGccqGGOaakcqWFbbqqdaqhaaWcbaGaemyzauMaemyzaugabaGaem4AaSgaaOGaeiykaKYaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqWFbbqqdaqhaaWcbaGaemyzauMaemiAaGgabaGaem4AaSgaaaGcbaGae8xqae0aa0baaSqaaiabdIgaOjabdkgaIbqaaiabdogaJbaakiabg2da9iab=feabnaaDaaaleaacqWGObaAcuWGIbGygaqbaaqaaiabdUgaRbaakiabgkHiTiab=feabnaaDaaaleaacqWGObaAcqWGLbqzaeaacqWGRbWAaaGccqGGOaakcqWFbbqqdaqhaaWcbaGaemyzauMaemyzaugabaGaem4AaSgaaOGaeiykaKYaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqWFbbqqdaqhaaWcbaGaemyzauMafmOyaiMbauaaaeaacqWGRbWAaaaakeaacqWFbbqqdaqhaaWcbaGafmOyaiMbauaacqWGObaAaeaacqWGJbWyaaGccqGH9aqpcqWFbbqqdaqhaaWcbaGafmOyaiMbauaacqWGObaAaeaacqWGRbWAaaGccqGHsislcqWFbbqqdaqhaaWcbaGafmOyaiMbauaacqWGLbqzaeaacqWGRbWAaaGccqGGOaakcqWFbbqqdaqhaaWcbaGaemyzauMaemyzaugabaGaem4AaSgaaOGaeiykaKYaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqWFbbqqdaqhaaWcbaGaemyzauMaemiAaGgabaGaem4AaSgaaaGcbaGae8xqae0aa0baaSqaaiqbdkgaIzaafaGafmOyaiMbauaaaeaacqWGJbWyaaGccqGH9aqpcqWFbbqqdaqhaaWcbaGafmOyaiMbauaacuWGIbGygaqbaaqaaiabdUgaRbaakiabgkHiTiab=feabnaaDaaaleaacuWGIbGygaqbaiabdwgaLbqaaiabdUgaRbaakiabcIcaOiab=feabnaaDaaaleaacqWGLbqzcqWGLbqzaeaacqWGRbWAaaGccqGGPaqkdaahaaWcbeqaaiabgkHiTiabigdaXaaakiab=feabnaaDaaaleaacqWGLbqzcuWGIbGygaqbaaqaaiabdUgaRbaaaaaaaa@AD20@

The set of equations Ak = 0 is somewhat of an oversimplification of the process of eliminating extraneous nodes. More precisely, the right hand side corresponding to each extraneous node is set to 0. To take a simple example, assume it is desired to eliminate a variable x2. If the row in the matrix Ak corresponding to x2 is equal to -3x1 + 4x2 - 2x3, then setting this row equal to 0 and solving for x2 yields x2 = 1/4(3x1 + 2x3). 1/4(3x1 + 2x3) is then substituted for x2 wherever it appears in any other equation, effectively eliminating x2 from the system of equations.

Total system matrix

The total system matrix comprises a combination of the preexisting FEM heart matrix (Ao), the FEM heart to torso coupling matrix Ac, and the FEM torso matrix At.

The FEM torso matrix At is based on a voxel based mesh of the torso. If the voxel torso mesh did not exist a priori, it may be generated by extending the preliminary mesh until it approaches the outer torso, as indicated by the green filled circles shown in Figure 3. The outer torso nodes may be found in a number of ways. In the approach adopted, the original torso surface nodes were multiply interpolated with two dimensional Legendre polynomials, resulting in a set of surface indicator functions analogous to the spherical harmonic surface indicator function described with reference to generation of the boundary surface.

The FEM stiffness matrix At for the voxel matrix was computed once again with a Galerkin type scheme with linear basis functions. For voxel meshes, this type of FEM matrix may be generated very quickly, on the order of a few seconds assuming reasonable computer processing speed and memory.

The FEM matrix (As) for the entire system is:

A s = [ A i i o A i h o A h i o A h h c + A h h o A h b c A b h c A b b c + A b b t A b g t A g b t A g g t ] ( 5 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaahaaWcbeqaaiabdohaZbaakiabg2da9maadmaabaqbaeqabqabaaaaaeaacqWFbbqqdaqhaaWcbaGaemyAaKMaemyAaKgabaGaem4Ba8gaaaGcbaGae8xqae0aa0baaSqaaiabdMgaPjabdIgaObqaaiabd+gaVbaaaOqaaaqaaaqaaiab=feabnaaDaaaleaacqWGObaAcqWGPbqAaeaacqWGVbWBaaaakeaacqWFbbqqdaqhaaWcbaGaemiAaGMaemiAaGgabaGaem4yamgaaOGaey4kaSIae8xqae0aa0baaSqaaiabdIgaOjabdIgaObqaaiabd+gaVbaaaOqaaiab=feabnaaDaaaleaacqWGObaAcuWGIbGygaqbaaqaaiabdogaJbaaaOqaaaqaaaqaaiab=feabnaaDaaaleaacuWGIbGygaqbaiabdIgaObqaaiabdogaJbaaaOqaaiab=feabnaaDaaaleaacuWGIbGygaqbaiqbdkgaIzaafaaabaGaem4yamgaaOGaey4kaSIae8xqae0aa0baaSqaaiqbdkgaIzaafaGafmOyaiMbauaaaeaacqWG0baDaaaakeaacqWFbbqqdaqhaaWcbaGafmOyaiMbauaacqWGNbWzaeaacqWG0baDaaaakeaaaeaaaeaacqWFbbqqdaqhaaWcbaGaem4zaCMafmOyaiMbauaaaeaacqWG0baDaaaakeaacqWFbbqqdaqhaaWcbaGaem4zaCMaem4zaCgabaGaemiDaqhaaaaaaOGaay5waiaaw2faaiaaxMaacaWLjaWaaeWaaeaacqaI1aqnaiaawIcacaGLPaaaaaa@77E0@

where the subscript g denotes the set of torso mesh nodes excluding the boundary nodes, and the subscript i denotes the set of original heart nodes excluding the outer surface heart nodes Xh. In this equation, the preexisting FEM heart matrix Ao has been partitioned into components: (i) that connect interior heart nodes to other interior heart nodes ( A i i o MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaqhaaWcbaGaemyAaKMaemyAaKgabaGaem4Ba8gaaaaa@3207@ ); (ii) that connect interior heart nodes to nodes on the heart surface ( A h i o MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaqhaaWcbaGaemiAaGMaemyAaKgabaGaem4Ba8gaaaaa@3205@ and A i h o MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaqhaaWcbaGaemyAaKMaemiAaGgabaGaem4Ba8gaaaaa@3205@ ); and (iii) that connect outer heart nodes to other outer heart nodes ( A h h o MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaqhaaWcbaGaemiAaGMaemiAaGgabaGaem4Ba8gaaaaa@3203@ ). Similarly, Ac and At have been partitioned into portions: (i) that connect "internal nodes" (i.e. nodes that are not a part of any other matrix) to other "internal nodes"; (ii) that connect "internal nodes" to "overlapping nodes" (i.e. nodes that appear in two matrices); and (iii) that that connect "overlapping nodes" to "overlapping nodes."

For example, A h h o MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaqhaaWcbaGaemiAaGMaemiAaGgabaGaem4Ba8gaaaaa@3203@ represents the coupling between the outer heart surface nodes, which are overlapping nodes, to themselves in the matrix Ao. Similarly, A h h c MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaqhaaWcbaGaemiAaGMaemiAaGgabaGaem4yamgaaaaa@31EB@ represents the coupling between the outer heart surface nodes to themselves in the matrix Ac. The matrix entry in As for the coupling of the outer heart nodes to themselves is thus the combined entry A h h o MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaqhaaWcbaGaemiAaGMaemiAaGgabaGaem4Ba8gaaaaa@3203@ + A h h c MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaieqacqWFbbqqdaqhaaWcbaGaemiAaGMaemiAaGgabaGaem4yamgaaaaa@31EB@ .

The discretized version of equation (1) is As V e = Ao, rhs V t, where V e and V t are the vectors corresponding to the continuous variables in equation (1), and Ao, rhs is the original, pre-computed stiffness matrix corresponding to the right hand side of equation (1). Ao, rhs and V t must be appropriately augmented with 0's to match the size of As V e. The boundary conditions at the heart/torso boundary, that Ve is continuous and the normal component of the current is continuous[3], are necessarily enforced by the finite element formulation used to create As.


The above described scheme was tested as follows. A spherical shell, analogous to the outer heart surface, was coupled to a surrounding spherical "torso". The potentials on the inner spherical shell were set equal to those which would have occurred had there been an ideal current dipole, normalized to unit strength, located at the sphere's origin and oriented along the z axis. Given this potential distribution, the potentials throughout the spherical "torso" were computed and compared with the analytic solution[8]:

V ( r , ϕ ) = cos ( ϕ ) ( 1 r 2 + 2 R r 3 ) ( 6 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGwbGvcqGGOaakcqWGYbGCcqGGSaaliiGacqWFvpGAcqGGPaqkcqGH9aqpcyGGJbWycqGGVbWBcqGGZbWCcqGGOaakcqWFvpGAieaacqGFPaqkcqGHxiIkdaqadaqaamaalaaabaGaeGymaedabaGaemOCai3aaWbaaSqabeaacqaIYaGmaaaaaOGaey4kaSYaaSaaaeaacqaIYaGmcqWGsbGuaeaacqWGYbGCdaahaaWcbeqaaiabiodaZaaaaaaakiaawIcacaGLPaaacaWLjaGaaCzcamaabmaabaGaeGOnaydacaGLOaGaayzkaaaaaa@4BD9@

where Φ is the polar angle (as before), R is the radius of the outer spherical torso surface, and r is the radius to an observation point. For the test problem, R was set at 50. The radius of the spherical shell/heart surface was set at 40, roughly the same size as the heart shown in Figure 1. To create the potential distribution equivalent to a dipole oriented along the z axis, the potentials on this sphere/heart surface were thus set at:

V ( 40 , ϕ ) = cos ( ϕ ) ( 1 40 2 + 2 50 40 3 ) . ( 7 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGwbGvcqGGOaakcqaI0aancqaIWaamcqGGSaaliiGacqWFvpGAcqGGPaqkcqGH9aqpcyGGJbWycqGGVbWBcqGGZbWCcqGGOaakcqWFvpGAcqGGPaqkcqGHxiIkdaqadaqaamaalaaabaGaeGymaedabaGaeGinaqJaeGimaaZaaWbaaSqabeaacqaIYaGmaaaaaOGaey4kaSYaaSaaaeaacqaIYaGmcqGHxiIkcqaI1aqncqaIWaamaeaacqaI0aancqaIWaamdaahaaWcbeqaaiabiodaZaaaaaaakiaawIcacaGLPaaacqGGUaGlcaWLjaGaaCzcamaabmaabaGaeG4naCdacaGLOaGaayzkaaaaaa@4FC7@

A voxel mesh, akin to the torso mesh shown in Figure 3, was generated to fill in the space between the inner and outer spherical surfaces, and the corresponding FEM matrix At was computed. Next, a coupling matrix Ac to link the torso to the nodes on the spherical shell/heart surface was created. The system of equations As *V c = 0 was constructed, where As is the total system FEM matrix and V c is the solution. The spherical shell/heart surface nodes were set equal to their corresponding potential V(40, Φ), resulting in an augmented system A'*V' = b, where V(i)' = 1 where i represents the spherical shell/heart surface nodes. This system was solved using the biconjugate gradient stable method, after preconditioning by scaling each row in the matrix so that its diagonal element was equal to 1.

Relative error was defined as:

R E = ( V c V a ) 2 V a 2 ( 8 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGsbGucqWGfbqrcqGH9aqpdaGcaaqaamaalaaabaWaaabqaeaacqGGOaakcqWGwbGvdaWgaaWcbaGaem4yamgabeaakiabgkHiTiabdAfawnaaBaaaleaacqWGHbqyaeqaaOGaeiykaKYaaWbaaSqabeaacqaIYaGmaaaabeqab0GaeyyeIuoaaOqaamaaqaeabaGaemOvay1aa0baaSqaaiabdggaHbqaaiabikdaYaaaaeqabeqdcqGHris5aaaaaSqabaGccaWLjaGaaCzcamaabmaabaGaeGioaGdacaGLOaGaayzkaaaaaa@44DB@

where Vc is the computed solution (excluding the spherical shell/heart nodes) and Va is the analytical solution. Both the number of nodes on the spherical shell/heart surface and the number of nodes in the "torso" were varied.

Heart beat

The heart in Figure 1, which is based on the Auckland canine model[2], was meshed[1] and oriented within a human shaped torso, as shown. A cellular automata like scheme was used to generate an activation and repolarization sequence throughout the heart, resulting in a corresponding sequence of transmembrane potentials V t . Maximum V t was 100 mV, with the resting potential set at 0 mV. Initial left ventricular and right ventricular activation sites were chosen roughly in accord with Selvester et al.[9] but were adjusted to generate better matches of epicardial potentials with recordings made by Spach et al.[10] in intact chimpanzees.

A propagation sequence was generated by a simple eikonal like scheme, according to which propagation velocity (v) was assumed to be proportional to the square root of bulk cardiac tissue conductivity[11]. The eikonal like scheme was implemented by calculating a propagation time between connected nodes in the heart mesh, selecting initial activation sites, and then stepping through time to activate the remaining heart nodes. The activation time of a heart node was set equal to the activation time of its nearest neighbor plus the approximate propagation time from the heart node to this nearest neighbor. The propagation time between neighboring nodes was set equal to v/d, where v is velocity (proportional to the square root of the conductance between nodes) and d is the distance between nodes.

Repolarization was initiated from the right ventricular breakthrough area[10] and allowed to spread like a wave across the epicardium in the manner in which activation was simulated. Repolarization was then allowed to proceed inward transmurally toward the epicardium based simply on relative transmural depth. In other words, repolarization proceeded from the epicardium inward. Although this oversimplified repolarization sequence is not necessarily strictly physiologically correct[12], it was chosen because of convenience and because it produced reasonable intra-epicardial gradients and epicardial/endocardial transmembrane gradients that resulted in realistic T waves.

Within the heart muscle: (i) intracellular conductivity along and transverse to fiber direction was set equal to 1 and 1/6, respectively; (ii). extracellular conductivity along and transverse to fiber direction was set equal to 1 and 1/3, respectively. Ventricular blood and torso conductivity were set at 5 and 2.5, respectively.

The heart mesh had approximately 270000 nodes and the torso mesh had approximately 120000 nodes. The system of equations was solved with Matlab's preconditioned conjugate gradient routine. The preconditioner was an incomplete Cholesky (0 fill in) factorization of the system matrix.


For the sphere test, Figure 7 shows relative error as a function of the number of total nodes, i.e. heart surface nodes (inner spherical shell) and torso nodes (between inner and outer spherical shells). The error drops rapidly to somewhat over 4% but convergence is markedly slower thereafter. It should be noted that different error curves were obtained depending on the relative balance between heart surface nodes and torso nodes; the curve shown in Figure 7 is representative of the general pattern of fast convergence with relatively few nodes and slow convergence after that. If the number of heart surface nodes (inner spherical shell) was held constant and the number of torso nodes increased, convergence was not necessarily monotonic.

Figure 7
figure 7

Results of the sphere test. Convergence is not monotonic and is relatively slow after approximately a 4% error has been reached.

Figure 1 shows the resulting simulated electrocardiograms at various locations roughly corresponding to the 6 precordial leads in the standard 12 lead electrocardiogram. The electrograms have a reasonably realistic shape. The peak to peak amplitude of approximately 4–5 mV is in line with the approximately 4.5 mV voltage difference across the torso in a body surface mapping study of chimpanzees[10]. The simulated maximum voltage difference across the epicardium (which occurred approximately 40 ms into the QRS complex) was approximately 13 mV, which is similar to the maximum epicardial potential difference reported by Spach et al. 28 ms into the QRS onset of the chimpanzee QRS complex.

The top panel in Figure 8 is a cross sectional view through the torso, along the plane indicated in the lower panel, which also shows the potentials on the torso surface and the outer heart surface. (In the lower panel, the colorbar is scaled to the torso potentials; the larger magnitude maximum and minimum heart potentials are saturated in this color scale and are simply red and blue, respectively.) The positive and negative tissue on the outer heart surface is faithfully projected onto the torso surface (bottom panel), as desired. The transverse slice shows a potential jump from the depolarized outer heart surface to the adjacent torso volume, as theory predicts[13]. The potential distribution across the torso slice appears to be reasonable.

Figure 8
figure 8

Computed potentials. Potentials on the heart and torso are shown in the lower panel and potentials through the cross section indicated by arrows are shown in the upper panel. In the lower panel, the heart potentials, which are greater than the torso potentials, saturate the color scheme.


This work describes a convenient scheme for coupling a pre-existing, possibly curvilinear, finite element mesh to a voxel mesh. The scheme is easy to code and generates a complete system matrix (i.e. the coupling and torso matrices along with final matrix assembly) fairly quickly, on the order of 10's of seconds for the heart/torso coupling example problem on a 2.80 GHz Xeon processor running Microsoft Windows XP with 3 GB of DRAM. The relatively good speed is helpful in cases, such as in the example problem, where some experimentation is required to fix the orientation of the preexisting and voxel meshes.

In this problem, the only user supervision involved orienting the heart within the torso. The rest of the process, namely the generation of coupling and torso matrices and assembly of the final matrix, was automatic.

The disadvantages of the scheme are that it appears to converge relatively slowly after a given degree of accuracy has been reached. It also may converge non-monotonically if only the number of nodes in the voxel mesh is increased. Given the weight of advantages and disadvantages, the scheme may be most useful for prototyping, where high accuracy is not required at the outset and experimentation is required to align the preexisting and voxel meshes.

One possible reason for the slow convergence is that the boundary surface Sb can not perfectly fit the boundary nodes because this surface can not follow the sharp corners between some of the adjacent boundary nodes. Indeed, for this reason, the error of the fit between Sb and the boundary nodes decreases slowly as the number of torso nodes is increased, which mirrors the overall slow convergence of the scheme.

There are a few possibilities for improving the scheme. First, instead of trying to fit the boundary surface to all of the boundary nodes (Figure 3), the corner boundary nodes (at the vertices of the rectangle defined blue red circles in Figure 3) could be excluded from the surface fit, yielding a better surface fit to the non-corner boundary nodes. These corner nodes could be coupled to the boundary surface through tetrahedral elements.

Another possible improvement would be to generate the coupling matrix with a very high spatial resolution mesh, and then couple this to a lower spatial resolution voxel mesh by removing extraneous nodes from the higher resolution mesh in the same manner that extraneous nodes were removed from the scaffolding matrix in the scheme described above. The coupling between higher and lower resolution meshes would be easy because both are voxel based meshes, so that it is easy to generate a low resolution mesh whose nodes are all within the set of the high resolution mesh nodes. For example, the high resolution mesh could be generated by creating nodes at the midpoints of the low resolution mesh, as shown in Figure 9. This mesh is non-conforming, due to the presence of nodes in the smaller elements that are not nodes of the neighboring larger elements. As an alternative to removing these nodes through elimination, which may cause numerical problems, this type of mesh may be treated in the manner described by Schimpf et al.[14]

Figure 9
figure 9

Alternative meshing scheme. A high resolution scaffolding mesh is coupled to a low resolution torso mesh.

Yet another alternative would be to couple the outer heart surface and the boundary surface with the boundary element method, which would avoid the requirement of eliminating extraneous nodes.

The present scheme may be applied to problems in bioengineering other than the heart/torso coupling problem. Further, the present scheme could be used in areas outside of bioengineering. However, extensions of the scheme to problems in which the preexisting mesh (e.g. the heart mesh) is not reasonably smooth may prove difficult because prismatic elements in the scaffolding matrix could be very irregularly shaped. Similarly, irregular node spacing on the preexisting mesh might also cause problems. If there are too few nodes on the outer surface of the preexisting mesh, then the scaffolding mesh might not have sufficient resolution for an accurate solution. If the physics of a problem require a very accurate representation of the geometry in the area of the boundary between the scaffolding mesh and voxel mesh, the present scheme will likely not prove suitable.

In summary, a scheme has been advanced that may be useful in situations where rapid computational speed and ease of use is a higher priority than great accuracy. Further, it may be possible to improve the convergence of the present scheme.


Ac – the desired coupling matrix between Xh and Xb'

Ak – the FEM scaffolding mesh matrix

Ao – the preexisting FEM heart matrix

As – the total system matrix, assembled from Ao, Ac, and At

At – the FEM torso matrix

Boundary nodes – the nodes within the preliminary voxel mesh that are closest to, but outside of, the outer heart surface

Boundary surface – a surface that is fitted to the boundary nodes

Outer heart surface – the outer surface of the preexisting heart mesh

Preliminary voxel mesh – a voxel mesh that is slightly larger than the heart mesh, which is embedded within the preliminary voxel mesh

Scaffolding mesh – a mesh that is constructed that fills in the volume between the outer heart surface and the boundary surface

Sh – outer heart surface

Sb – boundary surface

Xp – nodes in the preliminary voxel mesh

Xb – boundary nodes

Rh – radius of spherical harmonic surface that approximates Sh

X'p – Xp points translated to the origin of the heart coordinate system

Xh – the subset of Xp that lies within the heart

Xb' – nodes derived by projecting boundary nodes on to Sb

Xhp – the heart nodes as projected on the boundary surface; part of scaffolding mesh nodes

Xbp – the boundary nodes as projected on the heart surface; part of scaffolding mesh nodes


  1. Hopenfeld B: Spherical harmonic-based finite element meshing scheme for modeling current flow within the heart. Med Biol Eng Comput 2004, 42: 847–51. 10.1007/BF02345219

    Article  Google Scholar 

  2. Nielsen PM, LeGrice IJ, Smaill BH, Hunter PJ: Mathematical model of geometry and fibrous structure of the heart. Am J Physiol 1991,260(4 Pt2):H1365-H1378.

    Google Scholar 

  3. Fischer G, Tilg B, Modre R, Huiskamp GJ, Fetzer J, Rucker W, Wach PA: A bidomain model based BEM-FEM coupling formulation for anisotropic cardiac tissue. Ann Biomed Eng 2000,28(10):1229–43. 10.1114/1.1318927

    Article  Google Scholar 

  4. Thompson JF, Soni BK, Weatherill NP: Handbook of Grid Generation. Boca Raton, Florida:CRC; 1998.

    Google Scholar 

  5. Ingram DM, Causon DM, Mingham CG: Mathematics and Computers in Simulation. Mathematics and Computers in Simulation 2003, 61: 561–572. 10.1016/S0378-4754(02)00107-6

    Article  MathSciNet  Google Scholar 

  6. Wang Z, Hufford G, Przekwas : Adaptive Cartesian/adaptive prism grid generation for complex geometries. 35th Aerospace Sciences Meeting and Exhibit, AIAA: January 1997

  7. Roth BJ: Action potential propagation in a thick strand of cardiac muscle. Circ Res 1991, 68: 162–173.

    Article  Google Scholar 

  8. Wilson FN, Bayley RH: The electric field of an eccentric dipole in a homogeneous spherical conducting medium. Circulation 1950, 1: 84–92.

    Article  Google Scholar 

  9. Selvester RH, Wagner NB, Wagner GS: Ventricular excitation during percutaneous transluminal angioplasty of the left anterior descending coronary artery. Am J Cardiol 1988,62(16):1116–21. 10.1016/0002-9149(88)90560-7

    Article  Google Scholar 

  10. Spach MS, Barr RC, Lanning CF, Tucek PC: Origin of body surface QRS and T wave potentials from epicardial potential distributions in the intact chimpanzee. Circulation 1977,55(2):268–8.

    Article  Google Scholar 

  11. Keener JP: An eikonal-curvature equation for action potential propagation in myocardium. J Math Biology 1991, 29: 629–651. 10.1007/BF00163916

    Article  MathSciNet  Google Scholar 

  12. Yan GX, Shimizu W, Antzelevitch C: Characteristics of M Cells in Arterially Perfused Canine Left Ventricular Wedge Preparations. Circulation 1998, 98: 1921–1927.

    Article  Google Scholar 

  13. Colli Franzone P, Guerri L, Pennacchio M, Taccardi B: Spread of excitation in 3-d models of the anisotropic cardiac tissue. iii. effects of ventricular geometry and fiber structure on the potential distribution. Math Biosci 1998, 151: 51–98. 10.1016/S0025-5564(98)10004-4

    Article  MathSciNet  Google Scholar 

  14. Schimpf PH, Haynor DR, Kim Y: Object-free adaptive meshing in highly heterogeneous 3-D domains. Int J Biomed Comput 1996, 40: 209–25. 10.1016/0020-7101(95)01146-3

    Article  Google Scholar 

Download references

Author information

Authors and Affiliations


Corresponding author

Correspondence to Bruce Hopenfeld.

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Hopenfeld, B. A convenient scheme for coupling a finite element curvilinear mesh to a finite element voxel mesh: application to the heart. BioMed Eng OnLine 5, 60 (2006).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: