We study numerically the electromagnetic scattering properties of three dimensional (3D), arbitrary shaped composite dielectric metamaterials . Using integral equation techniques, we first derive a surface integral equation formulation which produces well -conditioned matrix equation. To solve the obtained integral equations, we apply a Galerkin scheme and choose the basis and testing functions as Rao-Wilton -Glisson defined on planar patches. We then develop an algorithm to speed up the matrix -vector multiplication s by employing the well -known method of moments (MoM) and the multi level fast mu ltipole algorithm on personal computer (PC) clusters. Some 3D numerical examples are presented to demonstrate the validity and accuracy of the proposed approach. INTEGRAL EQUATION FORMULATION Consider the case in figure 2 involving the problem of electromagnetic scattering by an arbitrary shaped, homogenous die lectric object B characterized by electrical properties ( , )ii which is surrounded by an infinite medium having ( , )oo . Using the equivalence principle, the inner and outer normal (N) and tangential (T) electric and magnetic fields can be considered as: N-EFIE(l)  ˆ ˆ ˆ J M Einc l l l l l l ln n n −  =−  T-EFIE(l)  tanˆ ˆ ˆ ˆ ˆ ˆ M J ( E )inc l l l l l l l l l ln n n n n n    −   =   N-MFIE(l)  1ˆ ˆ ˆ J M Hinc l l l l l l ln n n − −  =−  T-MFIE(l)  1 tanˆ ˆ ˆ ˆ ˆ ˆ J M ( H )inc l l l l l l l l l ln n n n n n −−   −   =   In (1) -(4), ˆlndenotes the inner unit normal for l=i and outer unit normal for l=o, l l l  = is the impedance of the host medium with   ( ), ( ) l i inner o outer= , (E , H )inc inc ll denote the incident electric and magnetic fields respectively, and Jl and Ml are equivalent electric and magnetic currents on the surface of the object which are relat ed to the total electric ( E) and magnetic ( H) fields as ˆ J ( ) H( )llr n r= ˆ M ( ) E( )llr r n= The operators and  operators a re defined as 2X [ X( ) .X( ) ]G ( , )l l lssjk r ds k ds r r r−       = +   X X( ) G ( , )lsr r r ds     =  Where l l lk = is the wavenumber and we define G(r,r ) ( )4ljk ReR r rR− = = − Which is the homogenous Green’s function . Among the various possibilities of combination of the normal and tangential electric and magnetic fields, the following mixed formulation is called the electric and magnetic current combined field integral equation (JMCFIE ) [1, 2], i.e.      22 12 11 22 2 11T- EFIE(l) ( 1) N- MFIE(l) T- MFIE(l) ( 1) N- EFIE(l)l l ll l l ll − == ==+−  −−  1 11 12 tan 21 22 tan[ ] [ ] Jˆ (E ) H [ ] [ ] Mˆ (H ) Einc inc N N N N ll inc inc N N N N llZZ n ZZ n −   −  =  +    For a composite dielectric object involving B -1 regions B 1, B2,…, B n with constitutive parameters ( , )nn located in a host medium 00( , ) , the surface integral equation can be written as 1 1 1 1 11, 12, tan 1, 0 0 0 1 1 1 21, 22, tan 2, 0 0 0ˆ [ ] [ ] ( (E ) H )J Mˆ [ ] [ ] ( (H ) E )B B B inc inc n N N n N N l l n B B B B B B inc inc n N N n N N l l n B B BZ Z n Z Z n − − − −  = = = − − −  = = =   −       =    +            DISCRETIZATION AND ITERATIVE SOLUTION In order to discretize the surface integral equation , one popular choice is the Rao -Wilton -Glisson (RWG) basis function [3], which discretizes the surface of the dielectric object into small elements each of which is triangle. Figure 3 depicts an RWG f unction on a pair of triangles. An RWG function An RWG function associated with the thi edge can be defined as ,2 f ( ) ,2 0,ii i i ii ii ilrT lr r T otherwise + + + − − −  = −     Where il is the length of the main edge, i is area of triangle iT associated with the edge, and i+ and i+ are respectively position vectors with respect to the free vertex of iT+ and iT− . Also note that in figure 3, t he total charg e is zero for each RWG function, ensuring that all edges of iT+ and iT− are free of line charges. Discretization using RWG functions leads to the following:  n n n n ,n ,n n n n n1TT + NK NT TK E J H MTK - NT TT + NKn n  −  =  In which discretized operators can be written as , , , ,1[ , , , ] ,, , , , ,G ( , ) G ( , )4TT p r q s p r q sp q r s n p q n p q p r n q s nS S S Sn p r q s p r q sjk l l jk l ldr r r dr dr r r dr− +−    =  −        , , , ,1[ , , , ] , , , , , , ,ˆˆ G ( , ) G ( , )42NT p r q s p r q sp q r s n p q n p q p r n q s p r nS S S Sn p r q s p r q sjk l l jk l ln dr r r dr n dr r r dr  − + − +    =   −          ,,[ , , , ] ,, ,,ˆ G ( , )4NK p r q sp q r s pq p r q s nSSn p r q slln dr r r dr++ =     ,,[ , , , ] ,, ,,G ( , )4TK p r q sp q r s pq p r q s nSSn p r q slldr r r dr++ =    In the far -field scattering electric field can be written as: 00 1 0 0 E J( ) M( )4jk r jk r sca SSjr e dr r e dr−  −  −     =− The discretization gives: 00 0 0J J 1 M M 0 0 11 2 J 0 , 11 2 M 0 , 11E ( ) ( )4 42 42NN jk r jk r sca n n n nSSnn N jk r n n n qSnq nq N jk r n n n qSnq nqjc v r e dr c v r e dr jlc e dr jk lc e dr  −  −  − == −+ == −+ ==    =−  =   −   NUMERICAL RESULTS The procedure described above is implemented for the solution composite dielectric metamaterials. In this section, the developed surface integral equation is analyzed by numerical examples. The process is implemented in MATLAB. Sphere with 1  radius (b) triangulation with  /6 mesh size We first consider a simple case which involves the problem of plane wave scattering by a uniform dielectric sphere in free space . The sphere has a radius of 02 , where 0 is the free space wavelength. NUMERICAL EXAMPLES Current density on sphere using numerical solution of 3J( , ) ( , ) E ( , )tscar t dt d r r r t t r t  −     = − − Readers can refer to the following paper s for in depth theoretical analysis and device simulations in MEMS and optics [ 4-14].