Introduction to the finite element method in electromagnetics download

Chapters 1 to 4 present general 2D and 3D static and dynamic formulations by the use of scalar and vector unknowns and adapted interpolations for the fields nodal, edge, face or volume.

Chapter 5 is dedicated to the presentation of different macroscopic behavior laws of materials and their implementation in a finite element context: anisotropy and hysteretic properties for magnetic sheets, iron losses, non-linear permanent magnets and superconductors. More specific formulations are then proposed: the modeling of thin regions when finite elements become misfit Chapter 6 , infinite domains by using geometrical transformations Chapter 7 , the coupling of 2D and 3D formulations with circuit equations Chapter 8 , taking into account the movement, particularly in the presence of Eddy currents Chapter 9 and an original approach for the treatment of geometrical symmetries when the sources are not symmetric Chapter Chapters 11 to 13 are devoted to coupled problems: magneto-thermal coupling for induction heating, magneto-mechanical coupling by introducing the notion of strong and weak coupling and magneto-hydrodynamical coupling focusing on electromagnetic instabilities in fluid conductors.

Chapter 14 presents different meshing methods in the context of electromagnetism presence of air and introduces self-adaptive mesh refinement procedures. Optimization techniques are then covered in Chapter 15, with the adaptation of deterministic and probabilistic methods to the numerical finite element environment.

Chapter 16 presents a variational approach of electromagnetism, showing how Maxwell equations are derived from thermodynamic principles. Download Introduction To The Finite Element Method In Electromagnetics books , This lecture is written primarily for the non-expert engineer or the undergraduate or graduate student who wants to learn, for the first time, the finite element method with applications to electromagnetics.

It is also designed for research engineers who have knowledge of other numerical techniques and want to familiarize themselves with the finite element method. Finite element method is a numerical method used to solve boundary-value problems characterized by a partial differential equation and a set of boundary conditions. Author Anastasis Polycarpou provides the reader with all information necessary to successfully apply the finite element method to one- and two-dimensional boundary-value problems in electromagnetics.

The book is accompanied by a number of codes written by the author in Matlab. These are the finite element codes that were used to generate most of the graphs presented in this book. Specifically, there are three Matlab codes for the one-dimensional case Chapter 1 and two Matlab codes for the two-dimensional case Chapter 2. The reader may execute these codes, modify certain parameters such as mesh size or object dimensions, and visualize the results.

Download Finite Elements Electromagnetics And Design books , Advanced topics of research in field computation are explored in this publication. Contributions have been sourced from international experts, ensuring a comprehensive specialist perspective. A unity of style has been achieved by the editor, who has specifically inserted appropriate cross-references throughout the volume, plus a single collected set of references at the end. The book provides a multi-faceted overview of the power and effectiveness of computation techniques in engineering electromagnetics.

In addition to examining recent and current developments, it is hoped that it will stimulate further research in the field. Download Frequency Domain Hybrid Finite Element Methods For Electromagnetics books , This book provides a brief overview of the popular Finite Element Method FEM and its hybrid versions for electromagnetics with applications to radar scattering, antennas and arrays, guided structures, microwave components, frequency selective surfaces, periodic media, and RF materials characterizations and related topics.

It starts by presenting concepts based on Hilbert and Sobolev spaces as well as Curl and Divergence spaces for generating matrices, useful in all engineering simulation methods.

It then proceeds to present applications of the finite element and finite element-boundary integral methods for scattering and radiation. Applications to periodic media, metamaterials and bandgap structures are also included.

The hybrid volume integral equation method for high contrast dielectrics and is presented for the first time. Another unique feature of the book is the inclusion of design optimization techniques and their integration within commercial numerical analysis packages for shape and material design.

To aid the reader with the method's utility, an entire chapter is devoted to two-dimensional problems. The book can be considered as an update on the latest developments since the publication of our earlier book Finite Element Method for Electromagnetics, IEEE Press, The latter is certainly complementary companion to this one.

In addition, the discretization of the domain allows the weak formulation of the problem to be applied to each element separately, thus allowing us to define distinct element values for material properties and sources. This offers generality and versatility to the method. Since the weak formulation of the problem, as it will be outlined in the following section, contains first-order derivatives of the primary unknown quantity, the chosen interpolation functions must be continuous within the element and at least once differentiable.

The use of a polynomial instead of any other type of function allows for the integrations, which are products of the weak formulation, to be evaluated more easily. Besides continuity and differentiability, another important requirement is that these polynomials must be complete. In other words, they must consist of all the lower order terms.

This is essential in order for the solution to be accurately represented by the interpolation functions inside an element. Let us now consider a finite element line segment , as illustrated in Figure 1. This element has coordinates x1e and x2e , which correspond to local nodes 1 and 2, respectively.

As illustrated in the figure, the master element has a fixed position along the natural coordinate axis. It is therefore easier for us to integrate a function on the natural coordinate system rather than on the regular coordinate system. In other words, by mapping an element onto the natural coordinate axis, the limits of integration involved in the weak formulation do not change every time a new element is considered. This would be the case if an integral were to be evaluated for elements residing on the regular coordinate axis.

Due to this observation, it is instructive that interpolation functions be derived based on the master element rather than the local element. In our case, there are two nodes per element and, therefore, it is necessary to use two interpolation functions. For the expression in 1. If this is not the case, then, the primary unknown quantity V will not be continuous across element boundaries.

This process is known as the method of weighted residual or weighted-residual method. Given the element shown in Figure 1. As mentioned in Section 1. The weighted residual is formed by moving all terms of the differential equation in 1.

To be able to solve for the unknowns, it is important to have as many independent weight functions as the number of nodes or unknowns. It is often preferred that the weight functions w x be identical to the interpolation or shape functions N x. From the weighted-integral equation in 1. The second element right-hand-side vector in 1. Remember that the entire domain is a collection of elements. The element matrices and vectors representing each finite element in the domain must be assembled according to the connectivity of the elements in order to form the global matrix system.

The assembly process will be explained in the following section. The entries of the element matrix K e are obtained by evaluating the integral in 1. However, in Section 1. Specifically, it can be shown from 1.

Solution, however, of the BVP mandates the assembly of all finite elements, based on the element connectivity information, and the formation of a global matrix system of linear equations. This global matrix system can be solved to obtain a numerical solution to the problem once the Dirichlet boundary conditions are imposed. The imposition of Dirichlet boundary conditions is the topic of the following section. To explain the process of assembly, consider Figure 1.

Each element e has two nodes with local node numbers 1 and 2. According to the weak formulation presented in the previous section and assuming linear finite elements, a set of two equations is produced for each element in the domain. Substituting this numerical solution into 1. Equations 1. The second global right-handside vector in 1.

On the contrary, V corresponds to the approximate numerical solution whose accuracy depends on the choice of subdomain interpolation functions. The first derivative of V is not necessarily continuous across elements, regardless of the fact that material parameters, such as dielectric constant, are the same for all elements.

The continuity of the secondary variable across elements, which is directly proportional to the first derivative of the primary variable, was never imposed in the finite element formulation presented in Section 1. In particular, the first derivative of V is constant over the length of a linear element and, as a result, there will be a step discontinuity in the distribution of the secondary variable.

This discontinuity will appear at element boundaries. However, as the number of elements in the domain is increased, this step discontinuity tends to decrease. The element connectivity information relates the local node numbers of an element to the corresponding global node numbers, as shown in Table 1. A global node number, starting from 1 ending at Nn where Nn is the total number of nodes in the finite element mesh , is uniquely assigned to each node.

Note that the global coefficient matrix K must be initialized to the e zero matrix before we begin the assembly process. Therefore, according to Table 1. To illustrate the assembly process, consider an example where we have only three elements in the domain, i. Each entry from the element right-hand-side vector is added to the corresponding entry of the global right-hand-side vector according to the element connectivity information tabulated in Table 1. Note that the global right-hand-side vector must be initialized to zero before the beginning of the assembly process.

Ignoring for now the contribution of element vector de and taking into account only the contribution by element vector f e , which is given by 1. In addition, for the filling of the global coefficient matrix, the Matlab command sparse was used in order to save on computer memory by avoiding storing the zero entries. Note that the majority of entries in the global coefficient matrix will be zero after the completion of the assembly process. We call such a matrix a sparse matrix.

The sparsity of the global coefficient matrix is attributed to the subdomain nature of the shape functions.

In other words, a given shape function, which corresponds to a specific node, is nonzero only inside the elements where this node belongs to. For the remaining elements in the discretized domain, this specific shape function is zero and, therefore, there is no interaction between the specific node and the associated nodes of those elements.

As a result, the corresponding entries of the global coefficient matrix will be zero. Prior to imposing boundary conditions, the global matrix system is singular and, thus, cannot be solved to obtain a unique solution.

A nonsingular matrix system is obtained after imposing the boundary conditions associated with a given BVP. The two types of boundary conditions that are discussed in this section are the Dirichlet or essential boundary conditions and the mixed boundary conditions. The Dirichlet boundary conditions involve only the primary unknown variable whereas the mixed boundary conditions involve both the primary unknown variable and its derivative.

Another type of boundary conditions is the Neumann boundary conditions which can be considered as a special case of the mixed boundary conditions. Thus, there is one degree of freedom dof or unknown per node. The matrix system in 1. The set of linear equations in 1. Similarly, after eliminating a given column, all the columns on the right must be shifted toward the left by one position. This method of imposing Dirichlet boundary conditions is known as the Method of Elimination because the algebraic equation that corresponds to the global node at which the Dirichlet boundary condition is imposed must be eliminated, thus reducing the size of the governing matrix system by one.

From a programming point of view, it is more convenient to number the global nodes of the finite element mesh in such a way that the nodes which correspond to Dirichlet boundary conditions appear last. Thus, using the method of elimination in imposing Dirichlet boundary conditions, it is not necessary every time we delete a row to have to shift the bottom rows upward. The same argument applies to columns as well. This will save computational time and make programming simpler and more straightforward.

For a Neumann boundary condition, i. As it will be shown in Chapter 2, a first-order absorbing boundary condition ABC , which is used in scattering and radiation problems to truncate the unbounded free space, is a form of a mixed boundary condition.

A detailed explanation on how a first-order ABC is implemented in a 2-D scattering problem will be provided in Chapter 2. The leftmost plate is maintained at a constant potential V0 whereas the rightmost plate is grounded. The analytical solution to the problem was presented in Section 1.

The exact analytical expressions for the electric potential and electric field that exist in the region between the two plates are given by 1. For a finite element simulation and comparison of the numerical solution with the exact analytical solution, it is required that certain parameters be defined. Thus, according to Section 1. All units are expressed in the metric system. Thus, based on the assembly process presented in Section 1.

Once this is done, the right-hand-side vector must be updated according to 1. To plot the electric potential at intermediate points requires the use of the interpolation or shape functions employed for each finite element. The electric field at any point inside an element is computed by taking the negative gradient of the electric potential given by 1. Applying 1. Unlike the electric potential, which is continuous across element boundaries, the electric field is discontinuous.

Use of higher order interpolation functions will result in a better representation of the electric field in the finite element domain. Using a four-element mesh, the electric potential over the entire domain is computed and plotted. The finite element solution is compared against the analytical solution obtained in Section 1.

The comparison is shown in Figure 1. As illustrated, the electric potential at the nodes of the finite element mesh matches perfectly the analytical solution, whereas at intermediate evaluation points there is a deviation between the two solutions.

The reason for this deviation stems from the fact that the numerical solution at intermediate points is an interpolation of the nodal values using linear shape functions. An acceptable representation of the numerical error between the finite element solution and the exact analytical solution is defined as the area bounded by the two curves, which are depicted in Figure 1.

From this figure, it is observed that by doubling the number of elements, which is equivalent to reducing the length of the elements to half, the percent error in the numerical solution, as compared to the exact analytical solution, is reduced by a factor of 4. This can be clearly seen from Table 1. TABLE 1. Comparing Figure 1. Note that by doubling the number of linear elements in the finite element domain, the error based on the L2 norm decreases by a factor of 4, which was also the case for the numerical percent error calculated based on the area bounded by the two solutions.

Either way of computing the numerical error is acceptable. Both methods provide a good indication of the accuracy of the finite element solution as compared to the exact analytical solution. Besides plotting the electric potential in the finite element domain, one may decide to plot the electric field as a function of the x-coordinate.

The exact analytical expression of the electric field is given by 1. The finite element solution, based on linear elements, is given by 1. As shown from this expression, the electric field is constant inside an element and certainly not necessarily continuous across element boundaries. Figure 1. Note that the electric field is a vector quantity and, for this problem, it has a direction along the positive x-axis.

As stated before, the electric field obtained from the numerical approach is shown to be constant over the element and discontinuous across element boundaries. As the finite element mesh becomes increasingly denser, the numerical solution approaches the exact analytical solution. In this section, we are introducing higher order interpolation functions, specifically quadratic and cubic, in order to more accurately represent the finite element solution within the discretized domain.

It is expected that the numerical error will be substantially reduced with the use of higher order elements as opposed to linear elements. A linear representation of the solution over an element requires the values of the primary unknown quantity at only two nodes, which coincide with the end nodes of the element. On the other hand, a quadratic representation of the solution over an element requires the values of the primary unknown quantity at three nodes instead of just two.

Two of these nodes coincide with the end nodes of the element whereas the third one must be an interior node. Although the third node could be chosen at any interior point, the most convenient choice is the midpoint of the element. These are also known as Lagrange shape functions. Substituting the respective quadratic shape functions into 1.

The element matrix and right-handside vector of the BVP at hand will be evaluated using quadratic elements in Section 1. Exercise 1. Two of these coincide with the end nodes of the element and the other two correspond to interior points. The element matrix and right-hand-side vector for the BVP at hand will be evaluated using cubic elements in Section 1. These will eventually be discarded after imposing the two Dirichlet boundary conditions at the end nodes of the domain.

To obtain the entries of the element coefficient matrix K e , the integral in 1. Doing this, the limits of integration remain always the same for all elements in the domain, i. Using quadratic elements, show that the element coefficient matrix is given by 1. Symbolic math packages, such as Maple, can be used to analytically obtain the entries of the element coefficient matrix and right-hand-side vector, thus avoiding unnecessary mathematical complexities by hand.

Follow the same formulation presented in Section 1. To avoid the evaluation of complicated integrals by hand, it is recommended that the symbolic math package called Maple be utilized. For a better representation of the unknown quantities over the computational domain, it is instructive that these quantities be evaluated at points other than the nodes of the mesh. For the specific BVP considered in this chapter, the primary unknown quantity is the electrostatic potential whereas the secondary unknown quantity is the electric field.

Therefore, the electrostatic potential can be easily plotted in terms of the x-coordinate by looping through all the elements one-by-one and evaluating its value, using 1.

It is interesting to observe that the numerical solution is identical to the exact analytical solution, and therefore, the two plots are indistinguishable. Thus, the numerical error is effectively zero since quadratic elements were used to interpolate a solution which is quadratic in nature [see Eq. This would not have been the case if the exact analytical solution were of higher order. The electric field is computed by taking the first derivative of the electrostatic potential in 1.

To evaluate the electric field in the domain, we must loop through all the elements and, for each element, evaluate the expression 1. Using ten evaluation points per element, the x-directed electric field between the two parallel plates is shown plotted in Figure 1. As was the case with the electrostatic potential, the finite element solution is identical to the exact analytical solution, thus resulting in zero numerical error.

The reason stems from the fact that the exact expression of the electric field within the two parallel plates is linear in nature [see Eq. It is worth emphasizing here that the numerical error is not identical to zero but close to the machine error of the computer used to perform the computations. The separation between the plates and the boundary conditions are maintained the same as before.

The finite element formulation for this type of problem is exactly the same as the one outlined for a uniform charge distribution with the only difference being the element right-handside vector.

Using the nonuniform charge distribution described by 1. The distribution of the governing electrostatic potential is plotted in Figures 1. In both figures, the exact analytical solution, given by 1. It can be clearly seen that the numerical solution closely approaches the exact solution as the number of quadratic elements in the mesh is increased.

The same conclusive remark applies to the electric field. However, the electric field is linearly interpolated because it is computed from the gradient of the electrostatic potential. In addition, there is a discontinuity of the electric field at element boundaries. This was the case with linear elements as well. This discontinuity of the secondary unknown quantity tends to become increasingly smaller as the number of elements increases.

Note that the underlined weak formulation and associated shape functions guarantee continuity of the primary unknown quantity across elements but not necessarily of the secondary unknown quantity. This is demonstrated by plotting the electric field distribution between the two parallel plates for a two-element mesh and a four-element mesh using quadratic shape functions. The corresponding graphs, together with the exact analytical solution, are illustrated in Figures 1.

Once again, the finite element solution, for both the electrostatic potential and the electric field, becomes more accurate as the number of quadratic elements increases. The accuracy of the numerical solution can be evaluated by computing the numerical error, based either on the area bounded between the two curves or the definition of L2 norm, as a function of the number of quadratic elements.

This was done here using the first definition of numerical error, i. This is referred to as the numerical percent error. A plot of the numerical percent error, as far as the computation of the electrostatic potential is concerned, as a function of the number of quadratic elements is depicted in Figure 1. Comparing this graph with the corresponding graph of Figure 1. The rate of convergence using quadratic shape functions can be more clearly seen from Table 1.

Show that 1. Remember that the electron charge distribution in the region between the two parallel plates has changed from a constant profile to a parabolic profile given by 1. Using quadratic elements and the formulation presented in Section 1. The only difference is the order of shape functions used to interpolate the primary unknown variable over an element. To avoid unnecessary repetitions, we will restrict our discussion on the development of the governing expansion equations describing the primary and secondary unknown variables over an element after the finite element solution has been obtained.

It will be left as an exercise for the reader to repeat, using cubic elements, the same procedure as was done in the previous section, using quadratic elements, for the computation of the numerical percent error between the exact and the finite element solutions and to compare the result with Figure 1. Solving the global matrix system after imposing the governing Dirichlet boundary conditions, a set of values for the primary unknown variable i.

These values correspond to the global nodes of the finite element mesh. For someone to plot the primary unknown variable with good enough resolution in the discretized domain, it is necessary that the primary unknown variable be evaluated at multiple points along each element. The cubic shape functions expressed in terms of x are obtained by substituting 1. As was done with quadratic elements, the electric field over a cubic element is computed by taking the first derivative of 1.

From this discussion, it becomes clear now that the electric field, which is the secondary unknown variable, is quadratic over a cubic element whereas the electrostatic potential, which is the primary unknown variable, is cubic.

Dirichlet boundary conditions are imposed on the two plates. The code computes both the electrostatic potential and the electric field within the domain of interest. The third code, FEM1Dqnucd, uses quadratic elements to solve the same problem but with a nonuniform, instead of uniform, charge distribution. All three codes are capable of computing the numerical error associated with the finite element solution. Certain parameters such as mesh size and separation of plates can be modified by the user.

Write a finite element code in Matlab to solve the same BVP considered in this chapter but using cubic elements instead of linear or quadratic elements. Assume that the electron charge distribution between the plates is given by 1. Compute and plot the numerical error, either as a percentage or using the L2 -norm definition, and compare with Tables 1.

Zienkiewicz, The Finite Element Method. New York: McGraw-Hill, Silvester and R. Ferrari, Finite Elements for Electrical Engineers, 2nd ed. London: Cambridge University Press, Volakis, A. Chatterjee, and L. Pelosi, R. Coccioli, and S. Boston: Artech House, Mikhlin, Variational Methods in Mathematical Physics. New York: Macmillan, Sadiku, Numerical Techniques in Electromagnetics, 2nd ed.

Cheng, Fundamentals of Engineering Electromagnetics. New York: Addison-Wesley, Edwards, Jr. Penney, Elementary Linear Algebra. New Jersey: PrenticeHall, Such problems usually involve a second-order differential equation of a single dependent variable that is subject to a set of boundary conditions.

These boundary conditions could be of the Dirichlet type, the Neumann type, or the mixed type. The domain of the problem is a 2-D geometry with an arbitrary shape. Thus, an accurate representation of the domain in the context of the FEM presumes discretization of the domain using the most appropriate shape of basic elements called the finite elements.

The most commonly used finite elements in two dimensions are the triangular and quadrilateral elements. Both types of elements will be used in this chapter to solve 2-D BVPs in electromagnetics.

Furthermore, higher order elements will be used to illustrate the improvement in accuracy of the numerical solution without necessarily increasing the number of elements in the finite element mesh. Comparing 2. Using the FEM, the first step is to accurately represent the physical domain of the problem by a set of basic shapes called the finite elements. The use of a rectangle, for example, as a basic finite element to discretize an irregular domain is certainly the simplest but not the most suitable choice because an assembly of rectangles cannot accurately represent the arbitrary geometrical shape of the domain.

In such a case, the discretization error is significant, as shown in Figure 2. This is illustrated graphically in Figure 2. The quadrilateral is another basic element that is commonly used in 2-D finite element analysis.

A coarse mesh of the irregular domain using quadrilateral elements is shown in Figure 2. As was the case with the triangular element, the quadrilateral element results in a smaller discretization error than the one caused by the use of the rectangular element.

Note that there are certain advantages in using triangular elements as compared to quadrilateral elements, and these advantages become increasingly important when using vector elements to solve electromagnetic problems.

Since the scope of this book is to help the reader understand the basics of the FEM and not to apply the method to advanced topics in electromagnetics, our Discretization error a b FIGURE 2. For vector elements and their application to electromagnetic problems, the reader is referred to the literature [1—5]. The numbering of the nodes directly affects the bandwidth of the global matrix.

As was indicated in Chapter 1, these interpolation functions must satisfy certain key requirements. First, they must guarantee continuity of the primary unknown quantity across interelement boundaries. The triangle consists of three vertices which correspond to the three nodes of the element.

A linear interpolation function spanning a triangle must be linear in two orthogonal directions. Thus, a triangle of arbitrary shape, such as the one shown in Figure 2. Each linear interpolation function corresponds to a triangle node. If the local nodes of the triangle are numbered in a counter-clockwise direction, then the Jacobian comes to be positive and, therefore, it can be derived directly from the area of the triangle.

It is important to emphasize here that these interpolation functions are not linearly independent. A similar notation applies to the other two subtriangles. For isoparametric elements, the same shape functions used to interpolate the primary unknown quantity inside an element are also used to interpolate the space coordinates x and y. The quadrilateral has four local nodes which are numbered in a counter-clockwise direction.

Knowing the solution at the four nodes of the element, the primary unknown quantity can be evaluated at any point inside the element by using the appropriate interpolation functions. In this section, we will construct bilinear interpolation functions for the master quadrilateral element. An isoparametric representation will be used to transform a function from the natural coordinate system to the xy-coordinate system and vice versa.

Show that the shape functions governing a bilinear quadrilateral element are given by 2. The weak formulation of this problem can be obtained by first constructing the weighted residual of 2. The element residual is formed by moving the right-hand side of 2. However, this is not the case, and therefore, the element residual r e is, in general, nonzero. Our objective is to minimize this element residual in a weighted sense. To achieve this, we must first multiply r e with a weight function w, then integrate the result over the area of the element, and finally, set the integral to zero.

The contour integral in 2. Comparing the first integral of 2. Substituting 2. This integral will be treated separately a bit later. Equation 2. For example, if the finite element mesh consists of linear triangular elements, this contour integral must be evaluated along the three edges of each triangle in a counter-clockwise direction.

However, it is important to realize that a nonboundary edge belongs to two neighboring triangles, as shown in Figure 2. As a result of this observation, evaluating the line integral in 2.

The reason for the opposite sign stems from the fact that the two outward unit vectors normal to the common edge of the two neighboring triangles, as shown in Figure 2. Notice that the line integrals in 2.

Consequently, the contribution of the line integral in 2. Thus, the only contribution of the line integral in 2. For interior edges, as was stated before, its contribution is zero.

In this section, we are going to consider two types of interpolation functions: one for the linear triangular element, and another for the bilinear quadrilateral element.

Higher order elements will be considered separately in Section 2. The determinant of the Jacobian matrix is equal to twice the area of the triangular element provided that the local node numbers of the triangle follow a counter-clockwise sense of numbering.

Thus, in forming the connectivity information array of the finite element mesh, it is instructive that the local nodes of each triangle be numbered in a counter-clockwise direction.

Using 2. In other words, instead of integrating over the triangular element on the regular coordinate system, it is more convenient that the integration be carried out on the master triangle which lies on the natural coordinate system.

Another element matrix that is part of the governing linear system of equations is matrix e T given by 2. Using the Jacobi transformation in 2. If, however, g is a function of the space coordinates x and y, then it has to be mapped to the natural coordinate system using 2. Here, it is assumed that g is constant inside the integral and, thus, the entries of the right-hand-side vector f e are given, according to 2. Note that 2.

The integral in 2. This is equivalent to having a right-hand-side vector pe equal to the zero vector 0. Now, for a generic mixed boundary condition, as given by 2. Consider the triangular element illustrated in Figure 2. Numbering the local nodes of the element as indicated in the figure, the integral in 2.

This is equivalent to subtracting the coefficients of u e1 , u e2 , and e e e u e3 from the matrix entries K 11 , K 12 , and K 13 , respectively, i. The element right-hand-side vector pe for a boundary element locally numbered as shown in Figure 2. Given a linear triangular element on the xy-plane whose nodes are numbered in a counter-clockwise direction, prove that the determinant of the Jacobian matrix is twice the area of the triangle.

Exercise 2. Using the Jacobi transformation, evaluate the double integral in 2. Follow the same procedure as was done e for M Repeat the exercise for matrix T e to show that its entries are given by 2. Using the chain rule of differentiation, the partial derivatives of the interpolation functions in 2. Using the Jacobi transformation introduced for triangular elements in 2. These expressions are given by 2. This integral must be evaluated numerically.

Although there are many numerical methods [8, 9] that can be used to evaluate a double integration, some of which include midpoint rule and trapezoid rule, the most appropriate and widely used approach in the FEM is Gauss quadrature. Table 2. For a complete table of Gauss weights and points up to and including a point quadrature, the reader is referred to [10]. It is also recommended that double precision arithmetic be used when implementing Gauss quadrature in a computer code. It is important, however, to emphasize here that Gauss weights and points for triangles are different from the ones tabulated in Table 2.

For more information on Gauss quadrature for triangular domains, refer to [10, 11].