Al-Khwarizmi Engineering Journal Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) Finite Element Based Solution of Laplace's Equation Applied to Electrical Activity of the Human Body Zainab T. Baqer Department of Electrical Engineering/ College of Engineering/ University of Baghdad Email: zainab_alissa@yahoo.com (Received 8 November 2009; Accepted 4 November 2010) Abstract Computer models are used in the study of electrocardiography to provide insight into physiological phenomena that are difficult to measure in the lab or in a clinical environment. The electrocardiogram is an important tool for the clinician in that it changes characteristically in a number of pathological conditions. Many illnesses can be detected by this measurement. By simulating the electrical activity of the heart one obtains a quantitative relationship between the electrocardiogram and different anomalies. Because of the inhomogeneous fibrous structure of the heart and the irregular geometries of the body, finite element method is used for studying the electrical properties of the heart. This work describes the implementation of the Conjugate Gradient iterative method for the solution of large linear equation systems resulting from the finite element method. A diagonal Jacobi preconditioner is used in order to accelerate the convergence. Gaussian elimination is also implemented and compared with the Precondition Conjugate Gradient (PCG) method and with the iterative method. Different types of matrix storage schemes are implemented such as the Compressed Sparse Row (CSR) to achieve better performance. In order to demonstrate the validity of the finite element analysis, the technique is adopted to solve Laplace's equation that describes the electrical activity of the human body with Dirichlet and Neumann boundary conditions. An automatic mesh generator is built using C++ programming language. Initially a complete finite element program is built to solve Laplace's equation. The same accuracy is obtained using these methods. The results show that the CSR format reduces computation time compared to the order format. The PCG method is better for the solution of large linear system (sparse matrices) than the Gaussian Elimination and back substitution method, while Gaussian elimination is better than iterative method. Keywords: Finite element, ECG, PCG, volume conductor and GE. 1. Introduction Bioelectric field problems can be found in a wide variety of biomedical applications which range from single cells, to organs, up to models which incorporate partial to full humane structure. In this work a class of direct and inverse volume conductor problems which arise in electrocardiography is studied. The solutions to these problems have applications to defibrillation studies, detection and location of arrhythmias, impedance imaging techniques, and localization and analysis of spontaneous brain activity (in case of electroencephalography) in epileptic patients. Furthermore, they can, in general, be used to estimate the electrical activity inside a volume conductor, either from potential measurements at an outer surface, or directly from the interior bioelectric sources [1]. 2. The Conduction System of the Heart The human heart basically consists of the left and right atria, and the left and right ventricles. Each of these parts includes myocardial tissue surrounding a cavity. The right atrium collects and stores deoxygenated blood from the body. The right atrium adjoins to the right ventricle through the tricuspid valve. The blood of the right ventricle is pumped through the pulmonary This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 38 arteries passing the pulmonary valve and finally reaches the lungs where it is oxygenated. The left atrium collects and stores the oxygenated blood from the lungs. The border between the left atrium and the left ventricle is the mitral valve. The blood of the left ventricle is pumped through the aorta into the body passing the aortic valves. The left and right ventricles are separated by the inter- ventricular septum, the left and right atria by the interatrial septum. Figure 1 illustrates the conduction system of the heart that controls these cardiac contractions. The figure shows: a) the sino-atrial (S-A) node in which the normal rhythmic self excitatory impulse is generated, b) the internal pathways that conduct the impulse from the S-A node to the atrio-ventricular (A- V) node, c) the A-V node in which the impulse from the atria is delayed before passing into the ventricles, d) the A-V bundle, which conducts the impulse from the atria into the ventricles, and e) the left and right bundles of Purkinje fibers, which conduct the cardiac impulse to all parts of the ventricles. Fig.1. The Conduction System of the Heart [2]. The SA node creates an impulse of electrical excitation that spreads across the right and left atria. This impulse initiates the depolarization of the nerves and muscles of both atria, causing the atria to contract and pump blood into the ventricles. Repolarization of the atria follows. The impulse then passes into the atrioventricular (A- V) node, which initiates the depolarization of the right and left ventricles, causing them to contract and force blood into the pulmonary and general circulations. The ventricle nerves and muscles then repolarize and the sequence begins again. Ion movement in heart muscle constitutes a current flow, which results in potential difference in the tissue outside the fibers and on the surface of the body .These potential differences could be measured by placing electrodes on the surface of the body (Fig.2) and then displaying the result as an ECG [3]. Fig.2. The Chest Leads [4]. 3. Bioelectric Volume Conductor A general volume conductor can be defined as a region of volume,  , which has conductivity, , and permittivity, , in which resides a source current , vI , where the ( v ) signifies per-unit volume. Solving a volume conductor problem means finding expressions for electrical field, E, the potential,  , everywhere within the volume,  , and/or on one of the bounded surfaces, i . The more general formulation in terms of the primary current sources within the heart described by Poisson’s equation for electrical conduction [1]: vI . in  …(1) One can define a surface bounding the region which includes the sources and recast the formulation in terms of information on that surface, yielding Laplac's equation because the distributions of voltages on the surface are solved instead of current sources within a volume [1]: This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 39 0.   in  …(2) with Neumann boundary condition: 0 on T …(3) and Dirichlet boundary condition: 0.  n on T ...(4) This particular formulation is known as the Cauchy problem for Laplace's equation. If the field is determined from the source and conductor (heart), the problem is called a direct problem. If the source (potential of the heart) is determined from the known field and conductor (potential of the part of the surface of the body), the problem is called an inverse problem [1]. 4. Main Program The flow chart of the main program can be shown in figure 3. The following sections describe the details of each part. 4.1. Mesh Generation For simplicity, the proposed model is defined on a simple two dimensional domain (fig.4) which shows a cross section of the heart. To find the potential distribution for this model, the region is divided into a number of finite elements (meshed) as illustrated in Fig.5. A mesh could be created either inside the program or generated elsewhere using third party software. In this work the two methods are adopted. 1) Automatic Mesh Generation is built from scratch. 2) Automatic Mesh Generation using ANSYS package 4.2. Automatic Mesh Generation This program performs a mesh generation of an arbitrary solution domain. A few points are given to determine the general configuration of the region. Then the program automatically generates triangular or quadrilateral elements. Triangular elements are chosen. The subroutine Input accepts the data which defines the solution region outline and the material zones. 0 .   bA Fig.3. The Flow Chart of the Main Program. This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 40 Fig.4. The Proposed Model (a Cross-Section of the Heart). . Fig.5. Mesh Configuration. The basic steps involved in a mesh generation are as follows:  Subdivide the solution region into few quadrilateral blocks,  Separately subdivide each block into elements,  Connect individual blocks Each step is explained as follows: (A) Definition of Blocks The solution region is subdivided into quadrilateral blocks (Fig.6). Sub-domains with different parameters must be represented by separate blocks. As input data, block topologies and the coordinates at eight points describing each block is specified. Each block is represented by an eight-node quadratic isoperimetric element (Fig.7), i.e. for example: Block number one, its eight nodes (1, 2, 3, 8, 7, 6, 5 and 4) Block number two, its eight nodes (5, 6, 7, 16, 15, 17, 9 and 8) Block number three, its eight nodes (9, 17, 15, 14, 13, 12, 11and 10) . . . and so on. With natural coordinate system ),(  , the x and y coordinates are represented as [6]:    8 1 ),(),( i ii xx  …(5) i i i yy ),(),( 8 1     …(6) Fig.6. Subdivision of the Solution Region Into Quadrilateral Blocks. Fig.7. Eight-Node Quadratic Block. where ),(  i is a shape function associated with node i, and ),( ii yx are the coordinates of node i defining the boundary of the quadrilateral block as shown in figure 8. The shape functions are expressed in terms of the quadrilateral or This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 41 parabolic isoparametric element shown in figure 9. They are given by [6]: for corner nodes, 7,5,3,1)1)(1)(1( 4 1  iiiiii  …(7) for midside nodes 8,6,4,2 )1)(2( 2 1 )1)(1( 2 1 2222   i iiiii  …(8) Fig.8. Typical Quadrilateral Block [6]. From the above equations any node can be created inside or on the boundary of each block using the shape function and the eight nodes defining the boundary of the quadrilateral block. The properties of the shape functions are [6]: 1. They satisfy the conditions 1),( 1    n i i …(9)       ji ji i ,0 ,1 ),(  …(10) 2. They became quadratic along element edges ).1,1(   The subroutine INPUT(Algorithm 1) reads the number of points defining the mesh NPOIN, the number of blocks NBLOCK, the element type NTYPE (triangular or quadrilateral element), the number of coordinate dimension NDIME (2-D or 3-D), the nodes defining each block (as described previously), and the coordinates of each node in the mesh. Usually these data are read only once and stored in files. It is impracticable to read these data whenever the program is run.      Fig.9. Eight-node Serendipity Element [6]. INPUT Algorithm 1.Read NP ,Number of coordinate points defining the solution region 2.Read NBLOCK, Number of blocks or zones 3.Read NTYPE, The type of elements into which the structure is to be subdivided 4.Read NDIM, The number of coordinate dimensions, NDIM=2 5.Do IELEM=1,. . ., NBLOCK 6.Do I=1,. . .,8 7.Read NL(IELEM,I),NL is the connectivity matrix 8.end Do (I) 9.end Do (IELEM) 10.Do J=1,. . ., NP 11.Do I=1,. . ., NDIM 12.Read COARD(J,I),coordinate of each node 13. end Do (I) 14. end Do (J) 15. Read NDIVX, NDIVY, the number of element subdivisions to be made in the  and directions, respectively for all the blocks 16.Read weighting factors ii WandW )()(  for all the blocks 17.end Algorithm 1 The INPUT Algorithm The following is an algorithm that estimates the shape function (using equations 7 and 8). This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 42 Evaluation of the shape function //************************************** //SHAPEF(  , ) //************************************** 1.SHAPEE[1]=0.25*(1- )*(1- )*(- - -1); 2.SHAPEE[2]=0.5*(1- * )*(1- ); 3.SHAPEE[3]=0.25*(1+ )*(1- )*( - -1); 4.SHAPEE[4]=0.5*(1- * )*(1+ ); 5.SHAPEE[5]=0.25*(1+ )*(1+ )*( + -1); 6.SHAPEE[6]=0.5*(1- * )*(1+ ); 7.SHAPEE[7]=0.25*(1- )*(1+ )*(- + -1); 8.SHAPEE[8]=0.5*(1- * )*(1- ); 9.End Algorithm 2 Evaluation of the Shape Function. (B) Subdivision of Each Block For each block, NDIVX and NDIVY, the number of element subdivisions to be made in the  and directions, respectively are specified. Also the weighting factors ii WandW )()(  are specified allowing for graded mesh within a block. As an example if block five is desired to be divided into three divisions in the x-axis direction with different weighting (fig 6), NDIVX have to be equal to three and the weighting factor as follows: W(1)=1, W(2)=1, W(3)=4. In specifying, NDIVX, NDIVY, ii WandW )()(  care must be taken to ensure that the subdivision along block interfaces (for adjacent blocks) is compatible.  and is initialized to -1 so that the natural coordinates are incremented according to [6]: FW W T i ii . )(2    …(11) FW W T i ii . )(2    …(12) where    NDIVX j j T WW 1 )(  …(13)    NDIVY j j T WW 1 )(  …(14)        elementsquadraticfor elementslinearfor F ,2 ,1 …(15) Three types of elements are permitted: a) Linear four-node quadrilateral elements, b) Linear three-node triangular elements, c) Quadratic eight-node isoparimetric elements. (C) Connection of individual Blocks After subdividing each block and numbering its nodal points separately, it is necessary to connect the blocks and have each node numbered uniquely. This is accomplished by comparing the coordinates of all nodal points and assigning the same number to all nodes having identical coordinates. That is the coordinates of node 1 is compared with all other nodes, and then node 2 with other nodes, etc., until all repeated nodes are eliminated. The basic blocks of the automatic mesh generator is illustrated in algorithm (3). The OUTPUT function prints out the coordinates of the nodes and the element topologies (Figs.10, 11). Automatic Mesh Generator 1.INPUT 2.Generate 2.1 Subdivide the blocks into quadrilateral element 2.2 Eliminate the repeated nodes at block interfaces 3. Triangle //divides each four-node quadrilateral element into two triangular elements. The subdivision is done across the shorter diagonal 4. OUTPUT// provide the coordinates of the nodes, element topologies and material property of the generated mesh. 5.end Algorithm 3 Main Blocks of the Automatic Mesh Generator. This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 43 Fig.10. Element Topologies. Fig.11. Coordinates of the Nodes. 4.3. ANSYS Mesh Generator The ANSYS package can be used for mesh generation (fig 5). The 6-node triangular element type is selected for the proposed model. The coordinates of each node and the connectivity matrix (the matrix that describe each element and its 6-nodes) could be taken. A procedure is implemented to compress the 6-nodes to 3-node triangular element to fit the model's requirement (algorithm 4). This procedure also creates a new connectivity matrix (now each element is described by only 3-node). Compression Algorithm //***************************************** //Renumber the 6-node triangular elements to 3- node elements //***************************************** 1.Construct a Matrix,"flag_node" to contain the new numbering of the nodes 2. N=0,the new number 3.Do IELEM =1,...,NELEM;loop for all elements 4.Do INODE =1,2,3 5.Check if flag_node not utilized do 6.Store the node number in flag_node; 7.Store the new number in the connectivity matrix 8.Save the Coordinates of the node 9.Increase the number 10.end of for INODE 11.end of for IELEM 12.end Algorithm 4 Compression Algorithm. 4.4. Calculation of the Element Coefficient Matrix The element coefficient matrix is constructed for each element using geometry information, element connectivity and nodal coordinates. The element eija of the coefficient matrix may be regarded as the coupling between nodes i and j for element e, xli , yli is the coordinate of node i its value can be calculated as shown in algorithm 5: Calculate the Element Coefficient Matrix Algorithm 1.Read from files the local coordinates for the three nodes of the element e (XL, YL) 2.P1[1]=YL[2]-YL[3],Construct P and Q vector 3.P1[2]=YL[3]-YL[1] 4.P1[3]=YL[1]-YL[2] 5.Q[1]=XL[3]-XL[2] 6.Q[2]=XL[1]-XL[3] 7.Q[3]=XL[2]-XL[1] //Compute the area of the element. 8.AREA=0.5*ABS(P1[2]*Q[3]-Q[2]*P1[3]); 9. Area QQPP a jijie ij *4   [6] 10.end Algorithm 5 Element Coefficient Matrix. This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 44 Methods of calculating the coefficient matrix can be subdivided into two broad categories [7]: ASSEMBLY and Element by Element. 1. Assembly Method In assembly methods, sparse global element matrices are built using the element coefficient matrices (algorithm 6). Element ija is the coupling between node i and j. It can be obtained using the fact that the potential distribution must be continuous across inter-element boundaries. The contribution to the i, j position in ][a comes from all elements containing nodes i and j. For example if elements 1 and 2 have node 1 in common then: 211 1 1111 aaa  For another example if node 4 belongs to element 1, 2 and 3; and it represent the second node for element 1 and third node for element 2 and element 3 hence, 333 2 33 1 2244 aaaa  and so on. Calculate the Global Element Coefficient Matrix Algorithm 1. Do I=1,..,NELEM;loop all elements 2.Read from files the local coordinates(XL, YL) for the three nodes of the element I 3.P[1]=YL[2]-YL[3],Construct P and Q vector 4.P[2]=YL[3]-YL[1] 5.P[3]=YL[1]-YL[2] 6.Q[1]=XL[3]-XL[2] 7.Q[2]=XL[1]-XL[3] 8.Q[3]=XL[2]-XL[1] //Compute the area of the element. 9.AREA=0.5*ABS(P1[2]*Q[3]-Q[2]*P1[3]); 10. Area QQPP a jijie ij *4   11.Save ij ea in a proper place in a global matrix according to the node number 12.end Algorithm 6 Global Element Coefficient Matrix Algorithm. 2. Element by Element Method In element by element methods [7], a global system matrix is never created (algorithm 7). Explicit or iterative solvers are used to solve the equations. There are many different iterative solvers and the choice of which one to use depends on the type of problem being solved. Element by element methods lead naturally to a parallel solution strategy that may be applied to all the general problem types. Whether to use assembly or to use element by element methods depends on performance. One should not be surprised if the preferred method changes from machine to machine, from problem to problem or from time to time as hardware characteristics change. Although some researchers disagree about which method is faster. Two methods are adopted in this work (assembly and element by element method). Element-by-Element Method 1.Do I=1,..,NELEMI; loop all elements 2.P1[1]=YL[2]-YL[3] 3.P1[2]=YL[3]-YL[1] 4.P1[3]=YL[1]-YL[2] 5.Q[1]=XL[3]-XL[2] 6.Q[2]=XL[1]-XL[3] 7.Q[3]=XL[2]-XL[1] 8.AREA=0.5*ABS(P1[2]*Q[3]-Q[2]*P1[3]) /Determine coefficient matrix for Element I 9.Do i=1..3 10.Do j=1..3 11.a[i][j]=(P1[i]*P1[j]+Q[i]*Q[j])/(4*AREA) 12.Find row and column number (IR,IC) of each node 13.If(a[i][j]!=0)Do 14.Store this value in proper place in a comp-matrix and store its row number in row matrix 15.end if 16.end for j 17.end for i 18.end for I 19.end Algorithm 7 Element-by-Element Method. 4.5. Matrix Storage Scheme Any finite element problem for which the element ‘stiffness’ matrices fit in cache can be computed in a fast and efficient manner. As soon as the cache memory size is exceeded, performance will quickly degrade and be dominated by the speed of the slower main memory. So, what is the largest problem size or the maximum number of finite elements Nelsc that can theoretically reside in cache memory? For an approximate calculation, only three quantities are required. These are the size of the cache Cs , the number of floating point values required for each finite element (FPel) and the number of bytes required to represent the floating point This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 45 precision (Bfpp). These are simply related in the expression: fppel s BFP C Nels *  …(16) To reduce storage requirements, various strategies have been developed. Typically, two dimensional problems can be fairly large before memory size becomes an issue. However with three dimensional problems, memory quickly becomes a problem "to perform analyses with thousands, if not millions of elements" and in such cases, it may be preferable to use element by element methods]. In this work the model is two- dimension but a general program is implemented. 4.6. Matrix Compression Algorithm For this purpose "reduce storage" and since element coefficient matrix A[N][N] is a symmetric positive definite sparse matrix Let Ki the number of non zero elements in the i-th column and let Kmax maximum(Ki ) i=0,…N-1 . In our implementation matrix A is stored using two Kmax* N matrices AC and ROWS where AC is a compressed version of the matrix A and ROWS contains the indices of the elements of A stored in AC (algorithm 8). In figure (12) an example demonstrates how this algorithm acts. As an example if the dimension of the element coefficient matrix is 10000*10000 and the maximum number of none zero elements in column is equal to 50 so instead of 10000*10000=100,000,000 storage element required 1000*50*2=100,000 storage elements are utilized. Matrix Compression Algorithm //Compress A[N][N] to AC[Kmax][N] 1.Do j=1,..,N 2. Do i=1,..,N 3.If(A[i][j]!=0)Do 4.Store row number, i in proper place in matrix ROW 5.Save A[i][j] in proper place in AC matrix 6.end of if statement in step 3 7.end for i loop 8.end for j loop 9.end of the algorithm Algorithm 8 Matrix Compression Algorithm.                                        00504 03203 20121 ]5][3[ , 00307 01202 10143 ]5][3[, 00300 00007 01002 10240 00103 ]5][5[ ROWS ACA Fig.12. The Compressed Version of Matrix A. 4.7. Equation Solution The finite element approximation of equation (2) can equivalently be expressed as a system of N equations with N unknowns N ...,,.........1 (the electrostatic potentials, for example). In matrix form, the above system can be written as: bAΦ where )( ijaA  is the global coefficient (stiffness) matrix For volume conductor problems, A contains all of the geometry and conductivity information of the model [1]. Researchers have different views for the conductivity. Some of them consider the conductivity varying with the dimension (x and y axis in the two dimension model) and with the type of the tissue [8]. As an example the conductivity of the bone differs from the conductivity of the fat and so on. Others for simplicity consider unit conductivity [9] and this conductivity is considered in this work. 4.8. Boundary Conditions For the finite element method, it turns out that the Neumann condition is very easy to apply; while the Dirichlet boundary condition ( 0 on T , equation 3) takes a little extra effort. The Nemann boundary condition ( 0.  n on T , equation 4) is satisfied automatically within the Galerkin and variational formulations. To apply the Dirichlet boundary condition directly, the )( ija matrix is modified and the ib vector by implementing the following steps such that the ith value of i is known (algorithm 9). This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 46 Applying-Dirichlet-Boundary Condition 1. Subtract from the ith member of the r.h.s. the product of ija and the known value of i (call it i _  ; this yields the new right hand side, . _   ijii abb 2. Zero the ith row and column of A: .0  jiij aa 3. Assign .1  iia 4. Set the jth member of the r.h.s. equal to i _  5. Continue for each Dirichlet condition. 6. Solve the augmented system, .   bA 7.end Algorithm 9 Applying - Dirichlet - Boundary Condition. For example, the following 3*3 system,                                   3 2 1 3 2 1 333231 232221 131211 b b b aaa aaa aaa where 3 is known, 3 . Following the algorithm described above, the following augmented system is obtained,                                      3 3232 3131 3 2 1 2221 1211 100 0 0 ab ab aa aa Solving this system yields the solution set, ),,( 321  . In the following sections different methods are implemented to solve this system. In the direct problem equation (2) is solved using the measured voltages on the surface of the heart to calculate the voltages at the surface of the torso. The inverse problems are formulated as using measurements on the surface of the torso and calculating the voltages on the surface of the heart. Three methods are used to solve these equations: Precondition Conjugate Gradient method (PCG), Gaussian Elimination (GE) and Back Substitution (BS) method and iteration method. In the following sections a brief description, implementation and results of each method would be shown. 4.9. Precondition Conjugate Gradient Method The PCG algorithm is used for solving systems of linear equations in the form Ax=b where A is an N*N (symmetric positive definite sparse matrix) [10]. In this work a diagonal PCG algorithm is considered. The preconditioner M is a vector simply obtained in the following way: M[i]=A[i][i]. The implementation of precondition conjugate gradient method is shown in algorithm 10. In the description of the algorithm, k defines the iteration count, r specifies the step length, and beta denotes the correction factor. In the implementation the norm of the residual, d is used for the convergence check. Precondition Conjugate Gradient Algorithm // Solve Ax=b 1.Choose x0,intial guess 2.r0 = b-Ax0 3. Solve Mz0 = r0 4.P0 = z0 5. d0 = vv_product(r0,r0) Do k=0,1,. . .Kmax 6.qk = mv_product(A,Pk) 7. alpha = vv_product(Pk,qk) 8. alpha = dk / alpha 9. xk+1 = xk + alpha Pk 10. rk+1 = rk – alpha qk 11. Solve Mzk+1 = rk+1 11. dk+1 = vv_product(zk+1,rk+1) 12. If (SQRT(dk+1) < Tolerance)Exist 13. beta = dk+1 / dk 14. Pk+1 = zk+1 + beta Pk 15 end Do 16 end the algorithm Algorithm 10 Precondition Conjugate Gradient The elapsed time for this algorithm is measured as shown in Fig. 13. Increasing the size of the matrix increases the consumed time. The number of iterations required to reach convergence is also increased with increasing the size of the matrix. It can be shown that the number of iterations to reach convergence is less than the size of the matrix. This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 47 4.10. Gaussian Elimination and Back Substitution Method The Gaussian elimination algorithm has three nested loops. Several versions of the algorithm exist, depending on the order in which the loops are arranged. Algorithm 11 shows one version of Gaussian elimination, which is implemented. It converts a system of linear equations Ax=b to a unit upper-triangular system Ux=y. The matrix U shares storage with A and overwrites the upper- triangular portion of A. The element A[k,j] computed on line 6 is actually U[k,j]. Similarly the element A[k,k] equal to 1 on line 8 is U[k,k]. This algorithm leads to the LU factorization of A as a product L*U. After the termination of the procedure, L is stored in the lower-triangular part of A, and U occupies the locations above the principle diagonal. Tim e ve rsus s ize of m atrix w ith e ps=.01 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 500 1000 1500 2000 2500 3000 Size of ma trix (N) T im e ( s e c ) Fig.13. Time Consumed by PCG Method With Different Size of a Matrix. For k varying from 0 to n-1, the Gaussian elimination procedure systematically eliminates variable x[k] from equation k+1 to n-1 so that the matrix of coefficients becomes upper-triangular. In the kth iteration of the outer loop, an appropriate multiple of the kth equation is subtracted from each of the equations k+1 to n-1. The multiples of the kth equation are selected such that the kth coefficient becomes zero in equations k+1 to n-1 eliminating x[k] from these equations. A typical computation of the Gaussian elimination procedure in the kth iteration of the outer loop is shown in Fig.14. The kth iteration does not involve any computation on rows 1 to k- 1 or columns 1 to k-1. Thus at this stage, only the lower-right (n-k)*(n-k) sub-matrix of A (the shaded portion in Fig.14) is computationally active. The sequential run time of the procedure for large n is approximately 3/2 3n [11]. Gaussian Elimination Algorithm //convert a matrix A to an upper triangular matrix *********************************** 1. begin 2. Do k=0,..,n-1 3. Do j = k+1,.., n-1 4. A[k,j] = A[k,j] / A[k,k] 5. y[k] = b[k] / A[k][k] 6. A[k,k]=1 7. Do i = k + 1,..,n-1 8. Do j = k + 1,.., n-1 9. A[i,j] = A[i,j] – A[i,k] *A[k,j]; 10. b[i] = b[i] –A[i,k] * y[k] 11. A[i,k]=0 12. end j loop 13. end i loop 14. end k loop 15.end Gaussian_Elimination Algorithm 11 Gaussian Elimination. After the full matrix A has been reduced to an upper-triangular matrix U with ones along the principal diagonal, a back-substitution is performed to determine the vector x. Algorithm 12 shows a sequential back-substitution for solving an upper-triangular system of equations U x=y . Fig.14. A Typical Computation in Gaussian Elimination. This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 48 Back Substitution Algorithm *********************************** 1.Do k = n-1,..,0 2.x[k] = y[k] 3.Do i = k-1,..,0 4.y[i] = y[i] – x[k]* U[i,k] 5.end i loop 6.end k loop 7.end Back_Sustitution Algorithm 12 Back Substitution. The consuming time by a Gaussian Elimination method to solve the system BXA  is measured with different sizes of the matrix A as shown in Fig.15. It is obvious that the elapsed time by this algorithm is very large compared to the PCG's time (Fig.13); so for the suggested model PCG method is more efficient than the GE method. Using GE or PCG may not affect in a number of problems. The constructed matrix in this work is positive definite symmetric sparse matrix makes the PCG method more suitable. Tim e versus size of the m atrix for GE m ethod 0 5 10 15 20 25 30 35 40 45 0 500 1000 1500 2000 2500 3000 size of the m atrix T im e ( s e c ) Fig.15. Time elapsed by the GE and Back Substitution Method. 4.11. Iteration Method This is the third method used in this work for the solution, at a node k in a mesh with n nodes [6] ik n ki,1i i kk k AV A 1 V    …(17) where node k is a free node, A is the element coefficient matrix and kV is the potential at node k. This equation can be applied iteratively to all the free nodes (where the potentials are unknown). The iteration process begins by setting the potentials of fixed nodes (where the potentials are prescribed or known) to their prescribed values and the potentials at the free nodes are equated to zero or to the potential average. The accuracy in this method is changed with the number of iterations and there is a noticeable difference between the results when the number of iterations is changed. The time consuming is large (Fig.16) compared to the previous methods (Fig.17) so this method is not practical for such problems. Tim e ve rsus m atr ix's size for ite ration m ethod 0 50 100 150 200 250 300 350 0 500 1000 1500 2000 2500 3000 s ize of the m atr ix T im e ( s e c ) Fig.16. Elapsed Time by the Iteration Method with Different Sizes for the Matrix. Tim e com parition betw e en GE, PCG and iteration m ethods 0 50 100 150 200 250 300 350 0 500 1000 1500 2000 2500 3000 size of the m atrix T im e ( s e c ) GE and Back substitution method PCG methd Iteration method Fig.17. A comparison between the Elapsed Time of PCG, GE and Iteration Methods. 4.12. Direct and Reverse Solution In a direct problem voltages are measured on the surface of the heart and used to calculate the voltages at the surface of the torso, as well as within the volume conductor of the thorax. This is This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 49 implemented in this work by dividing the ECG cycle into steps for every 20 msec. In the direct solution a measurement is obtained for the entire surface (fig.18). There are many issues that affect these results. In the implementation a homogenous tissue is considered, a cross-section (2-D) of the heart is proposed. A simple model is suggested to approximate the electrical activity of the heart. In the future the model will be improved (homogenous and 3D) so as to get more accurate results. The inverse problems are formulated as using measurements on the surface of the torso and calculating the voltages on the surface of the heart. Better results are obtained by increasing the number of measurements on the surface of the torso. Fig.18. Potential Distribution in the Human Heart During the Direct Phase (Numbers Represent the Potential in mV). 5. Conclusions In studying the electrical activity of the heart, a volume conductor approach is used to approximate the voltages on the surface of the heart using the voltages at the part of the surface of the torso. The results show that increasing the size of the surface (known values) will increase the accuracy of the results. The difference between the normal results and the measured results is due to the assumption of homogenous tissue and the simple 2-D model suggested approximating the thorax. There are hopes the studies in this field will continue to get more accurate results. The results will be improved if a 3D and homogenous model are suggested since this model is more close to the fact. An automatic mesh generator is built from scratch. By decreasing the size of the elements the result are more accurate. A comparison between different methods is performed. The results show that the PCG method is preferable for a fast system than the GE, BS and iteration methods, while the iteration method is the most time consuming method. In addition compression of sparse matrices increases the speed of the system. Notations AV Atrio-Ventricular CG Conjugate Gradient CSR Compressed Sparse Row COARD Coordinate E Electrical Field ECG Electrocardiography EEG Electroencephalography FEM Finite Element Method GE Gaussian Elimination ID Identification vI Current per unit volume LU Lower Upper L Lower U Upper Matr Matrix mv_ prod Matrix-Vector Product N Number of nodes NDIVX Number of Division in the-  Axis Direction NDIVY Number of Division in the  - Axis Direction NBLOCK Number of Blocks NPP Number of Prescribed (fixed) nodes NELEM Number of Elements PCG Precondition Conjugate Gradient PVM Parallel Virtual Machine SA Sino-Atrial V Potential vv_prod Vector-Vector Product Vk Potential at subspace k Vect Vector Wx(i) Weight of Division i in the  - Axis Direction This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ Zainab T. Baqer Al-Khwarizmi Engineering Journal, Vol. 6, No. 4, PP 37 - 51 (2010) 50 Wy(i) Weight of Division i in the  - Axis Direction 1-D One-Dimension 2-D Two-Dimension Greek letters  Volume  Conductivity  Permittivity  Potential e Potential at element e  Test Function  Surface i Basis function at node i i Interpolation function  , Intrinsic coordinates 6. References [1] Christopher R. Johnson."Adaptive Finite Element and Local Regularization Methods for the Inverse ECG Problem". Biology Society 17th Annual International Conference, pages 209–210. IEEE Press. Center for Scientific Computing and Imaging, Department of Computer Science, University of Utah, Salt Lake City, Utah 84112 US, 1995. [2] Julian Johnson, "Wireless Electrocardiogram" [3] Wael Namat Al_Sahar, "The design and implementation of a PC based ECG recording and interpretation system", M.Sc Thesis, University of Technology, Electronic Engineering Dept, 1999. [4] Pilkington TC, Plonsey R: "Engineering Contributions to Biophysical Electrocardiography", IEEE Press, John Wiley, New York, 1982. [5] Tirupathi R. Chandrupatla and Ashok D. Belegundu "Introduction to Finite Elements in Engineering", copyright ©1991 by Prentice-Halls Inc. [6] [Matthew N. O. Sadiku "Numerical Techniques in Electromagnetics"; 2nd edition, © 2001 by CRC Press LLC. [7] C.O.Moretti, J.B. Cavalcante Neto, T.N. Bittencourt and L.F. Martha,"A Parallel Environment for Three-Dimensional Finite Element Method Analysis", 2000. [8] O. Skipa, D. Farina, C. Kaltwasser, O.D¨ossel, W. R. Bauer "Fast Interaction ECG Simulation and Optimization-Based Reconstruction of Depolarization" Institute for Biomedical Technique, University of Karlsruhe (TH), Germany Institute for Biophysics, University at Wurzburg, Deutschland Oleg.Skipa@ibt.uni- karlsruhe.de, 2004. [9] M.Okajima, K.Doniwa, T. Yamana, K. Ohta, N.Suzumura, A. Iwata and K. Kamiya. "Individualized Modeling of the Heart and Torso for Finite Element Method to be used in Forward Calculation of Body Surface Map", IEEE, University, Toyoake, Nagoya Area 470-11, Japan, 1989. [10] Antonio d' Acierno, Giuseppe De Pietro, Antonio Giordano, "A Parallel Implementation of the Conjugate Gradient Method on the Meiko CS-2", IRSIP _ CNR, Napoli (Italy), 1996. [11] Ananth Grama, Anshul Gupta, George Karypis and Vipin Kumar "Introduction to Parallel Computing" 2nd edition published 2003. This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/ )2010( 51- 37، صفحة 4، العدد 6مجلة الخوارزمي الھندسیة المجلد زینب توفیق باقر 51 العنصر المحدد لحل معادلة البالس مطبقة على النشاط الكھربائي للجسم األنساني زینب توفیق باقر جامعة بغداد \كلیة الھندسة \قسم ھندسة الكھرباء zainab_alissa@yahoo.com: البرید األلكتروني الخالصة .بیئة سریریة ي لقَیْاس في المختبِر َأو فاإلى الظواھِر الفسلجیِة التي صعبة برؤیةلَتزوید تخدم لدراسة التخطیط الكھربائي لعمل القلب لنماذج الحاسبات تس ف العدید من األ. َتتغّیُر على نحو ممیز في عدد ِمْن الشروِط الباثولوجیِة انھّافي التخطیط الكھربائي لعمل القلب أداة مھمة للطبیب ن ان تكتش مراض ممك ین یح ِبللقلالنشاِط الكھربائِي بمحاكاة. بھذا المقیاس ة ب ة كّمی ائي و صُل على عالق یط الكھرب یاء التخط ویة أش ر س ة غی انس .مختلف ر متج ب الغی بب التركی بس .للقلب طریقة العنصر المحدد استخدمت لدراسة الخواص الكھربائیة, اللیفي للقلب والجسم الھندسي الغیر منتظم دد یصف ھذا العمل بناء طریقة التدرج المرافق ر المح ة العنص ن طریق ة م رة الناتج ة الكبی ادالت الخطی ة المع ل مجموع ة لح رط . التكراری ة الش ان طریق ت و . استخدمت لتسریع االقتراب) PCG(االولي لجاكوبي القطري ا بنی ة طریقة الحذف الكاوسي ایض ة التكراری ع الطریق ا وم ت معھ تخدام . قورن م اس د ت لق ر . انواع مختلفة من طرق خزن المصفوفات مثل طریقة الصف المتفرق المضغوط للحصول على افضل اداء ل العنص ابات لتحلی لغرض التحقق من دقة الحس ب االنس ,المحدد ائي لقل اط الكھرب ف النش ي تص الس الت ة الب روط اسلوب تقني تبني لحل معادل من ش ان (ان ض ان ) العالم لت ونیوم ي . درش ي بن بكة ال د ش مول ي . نفس الدقة حصل علیھا من ھذه الطرق. في البدایة برنامج العنصر المحدد الكامل بني لحل معادلة البالس. ++Cبأستخدام لغة برمجة رط االول ة الش طریق .ة الحذف الكاوسي و طریقة الحذف الكاوسي افضل من الطریقة التكراریةطریقھي االفضل في النظام الخطي الكبیر من PCGلجاكوبي القطري This page was created using Nitro PDF trial software. To purchase, go to http://www.nitropdf.com/ http://www.nitropdf.com/