Home

The µ-diff toolbox: user guide - mu-diff

1. OO0OO0OO0O0O OOO00 l N j l Ni A l l WOOO OOOO0 D0000 GOO00 1 1 1 1 o 5 10 15 20 o 50 100 15 0 x Angle of reception degree NE 8 at a 250 300 350 a Obstacles b Radar cross section Absolute value of the total field Real part of the total field c Absolute value of the total field d Real part of the total field Imaginary part of the total field 20 15 of 10 r 0 05 gt 0 o 5 at 0 05 15 0 1 20 20 15 10 5 o 5 10 15 20 x e Imaginary part of the total field Figure 4 5 Results obtained for a complex geometry of penetrable obstacles with k 1 and k 2k solved using a single layer potential representation of the field The point source is located on 0 0 in the center of the geometry 81 82 Appendix A List of the y diff functions alphabetical order e BenchmarkCalderon Examples Benchmark Example of solution of penetrable scattering using Calderon projectors e BenchmarkDirichlet Examples Benchmark Example of solution of sound soft scattering using different integral equations e BenchmarkNeumann Ezamples Benchmark Example of solution of sound hard scattering using different integral equations e BenchmarkPenetrable Examples Benchmark Example of solution of penetrable scattering using single layer potential BlockCalderonP ro jector IntOperators Dense Interface Block Calderon projector block IntOperators Dense Int
2. LineWidth LINEWIDTH set the line width to LINEWIDTH same as the plot function e PlotCirclesOnFigure zdata ZDATA set the zdata of the figure to ZDATA same as the plot function When drawing disks on a figure containing some values do not forget to set the zdata to the max value An example is given below figure 1 surf X Y U PlotCirclesOnFigure O a 1 zdata max max U If X and Y are obtained by meshgrid let us define the mask matrix S by 0 if X i 7 Y i j is outside the obstacles S j 4p if XG j Y i j Qp This class of matrices is handy when plotting the potentials and is actually used by the two func tions and The function MaskMatrixObstacles which provides the matrix M has the following syntax S MaskMatrixObstacles X Y O a Let us explain how to extract the boundary points In the same way as for the MaskMatrixObstacles function JBoundaryOfObstacles leads to a matrix G where if X i j Y i j is outside the obstacles G i j 40 if X i j Y i j Qp as Similarly G is computed by G BoundaryOfObstacles X Y O a 61 3 6 Examples available in y diff Some examples are available in the Examples directory and especially in Examples Benchmark where each file is independent and solves a different problem Dirichlet Neumann or penetrable They should also be launched to check that p dif
3. L Opulr x A ulr x 15 We also have Proposition 1 3 The trace and the normal derivative trace of the operators Z and M are given by the following relations A 4 1 1 Lp Lp Yo XA 57 M A 1 5 4 1 4 VW Lp 51 N p y MA DA where I is the identity operator and for x T p H T and X MAT the four boundary integral operators are defined by bi HPD gt HPD Lex fG M HPT gt HPD MAK day Geey AWAP LY N HPD HPD No f aaayyy M p x D HYD HAT DAS a Any GC YAYAT 1 6 Al along the user guide the boundary integral operators are written with a roman letter e g L whereas the volume integral operators are denoted by a calligraphic letter e g Z According to 21 Theorems 3 4 1 and 3 4 2 the boundary integral operators L and D are invertible providing k is not an irregular frequency More precisely we have Theorem 1 4 Let Fp Q resp Fy Q7 be the countable set of positive wavenumbers k accumulating at infinity such that the interior homogeneous Dirichlet resp Neumann problem 1 7 v 0 resp Onv 0 onT Da k v in Q7 admits non trivial solutions Then the operator L resp D realizes an isomorphism from H T into H T resp from H T into H7 T if and only if k Fp Q resp k Fy Q7 These irregular frequencies k of Fp Q resp of Fy Q7 are exactly the square roots
4. and for q 1 M with q p byq as the vector between the centers Og and Op gt bpa O40p bpa bpaql l Qpq Angle Ox1 bpq 29 O Figure 2 1 Illustration of the notations for two disks Q and Qg and a point x QF Furthermore any point x is described by its global polar coordinates r x Ox r x lex 0 x Angle Ox1 r x or by its polar coordinates in the orthonormal system associated with the obstacle Q5 with p l M rp x Opx r x ltp II p x Angle Ope1 rp x Let us now build a basis of L T to approximate the integral operators To this end we first construct a basis of L I associated with Q for p 1 M If the circle T has a radius one and is centered at the origin then a suitable basis of L is the spectral Fourier basis of functions e mez We adapt this basis to the general case where ap 4 1 by introducing on one hand the functions Y m mez defined on R by Vm Z Vx R Ym x e2 and on the other hand the functions y 1 lt p lt m mez given by Pmitp x _ ember Vp 1 M Vm Z Vx Ip P x 0 MEZ VET ai See For p 1 M the family Yfa mez forms an orthonormal basis of L Tp for the hermitian inner product p2 p VIERT hadseery f dab Jar p x To build a basis of L T we introduce the functions of L T as the union of these M families 0 if Vo q 1 M Vme Z olny dFP ye
5. a D21 anM2 p 7 by I 1 a N 2 anL 2 1 a D22 an G Maa Az b2 with l l bp 1 a u p an nu r For M obstacles the same procedure can be applied to obtain I where e the matrix A of size M x M is such that 1 NPA anl ifq lt M APY a EO E 1 35 1 a DP anM ifq gt Mp e the unknown vector y Yq 1 lt q lt m collects the contributions of the elementary densities _ Par ifl lt q lt Mp ea Ge Mpeg N e the right hand side b bp 1 lt p lt m is given by bp 1 aju r andnu p 1 4 The penetrable case 1 4 1 The boundary value problem Let assume that the obstacles are now penetrable with a contrast function n such that Is np ifxeEQ 0 otherwise 24 where np are constant for each obstacle The boundary value problem to solve is given by find the total field ur such that A k ur 0 int A k n ur 0 mn ur 0 onT Onur 0 onl ur ui is outgoing to Q7 1 4 2 An example of integral equation for the penetrable case To rewrite this problem as a boundary integral equation let the field ur be rewritten as ut u in OF ur i u in Q7 where ut and u are solutions of the coupled problems A k jut 0 in QF A k n u 0 in Q7 ut u r w p onT nut nu Ir Oyul p oI u is outgoing to Q7 Following the EFIE procedure let the two fields ut and u be expressed as a single layer
6. a 2 al Ar i 6 4 2 0 2 4 6 0 50 100 150 200 250 300 350 x Angle of reception degree a Obstacles b Radar cross section History of convergence of the GMRES k 1 without restart Precond Absolute value of the scattered field 2l al prn GMRES residual log 12 L L L L 1 1 f L j 1 2 3 4 5 6 a 8 9 Iteration number c History of convergence of the GMRES Absolute value of the total field 10 0 4 0 35 0 3 0 25 0 2 015 01 X 0 05 10 o 8 6 4 2 o 2 4 6 8 10 x A e Absolute value of the total field 0 06 0 05 0 04 0 03 l 0 02 b 0 01 10 8 4 4 2 o 2 4 6 8 10 a x d Absolute value of the scattered field Real part of the total field 19 0 35 03 0 25 02 015 0 1 0 05 E o 0 05 0 1 Sa 0 15 o 8 6 4 2 o 2 4 6 8 10 x f Real part of the total field Figure 4 4 Neumann boundary value problem with a point source emitter solved using the single scattering preconditioned integral equation In the respective order is shown the obstacles the radar cross section the history of convergence of the GMRES and then the absolute value of the scattered field the absolute and real part of the total field 75 4 3 Mixing Dirichlet and Neumann boundary conditions Let us consider the following situation where we mix Dirichlet and Neumann boundary condi tions The scatterer is composed of Mp sound soft and My M Mp sound hard obstacles leading
7. then frourierTruncat ion uses formula 2 16 with k p as the wavenumber The following options are moreover available e M modes FourierTruncation Min MIN To force a minimal value M_modes p is then either the min value between MIN and formula e M modes FourierTruncation Tol TOL The tolerance set by default to 10719 is then set to TOL 3 3 3 Incident waves Generalities Two different incident waves are available in the diff toolbox the plane wave and the point source wave The user can build his own incident wave They are all located in the directory PreProcessing IncidentWave As explained in section 2 1 4 a right hand side b is decomposed by blocks each of these blocks representing one obstacle b bp p 1 a A different condition can be applied on two different obstacles e g Dirichlet on Q and Neumann on N2 or a different integral equation can be considered on each obstacle e g EFIE on Q and MFIE on 2 To this end diff builds each block separately thanks to the function which computes the vector bp According to the input data the function builds one of the available right hand sides described in Table On the other hand the common function Incident Wave computes the whole vector b For all obstacles p the function BlockIncidentWavel is called and the whole vector is assembled In addition for all the incident waves an interface function is available These easy to use in
8. 0 a ones l size O 2 O a RemoveDisk O a X 0 Y 0 N_scat size 0O 2 4 4 3 Writing and solving the BIE using diff Let us solve the problem using a dense storage and first write the matrix of the system which consists in assembling the different operator into a larger matrix Splus SingleLayer 0O a M_modes k Sminus IntegralOperator O a M modes k_int 2x eye N_scat N_scat Nplus DnSingleLayer O a M modes k Nminus IntegralOperator O a M modes k_int 4xeye N_scat N_scat Identity eye size Nplus 78 Boundary integral operator A Splus Sminus 0 5 Identity Nplus 0 5 Identity Nminus cleaning memory clear Splus Sminus Nplus Nminus Identity mu_matrix Second we need to introduce the right hand side here a point source Uinc PointSource O a M_modes k XS DnUinc DnPointSource O a M_modes k XS B Uinc DnUinc clear Uinc DnUinc Finally we solve the system here using a direct solver The two densities pt and p are here also extracted to simplify rho A B 6 Extracting the exterior and interior densities rho_plus rho 1l sum_modes rho_minus rho sum_modes 1l end clear rho 4 4 4 Post processing Radar cross section We compute now the radar cross section which is based on the exterior single layer potential theta_RCS 0 360 theta_RCS_rad theta_RCS
9. 1 M_modes M_modesD M_modesN sRight hand side Uinc DnPlaneWave O a M_modes k beta_inc DnUinc DnPlaneWave O a M_modes k beta_inc B alphaxetaxUinc l alpha DnUinc 6 Assembling Assembling zeros 2 N_scat N_scat 76 Weight zeros 2 N_scat N_scat for p 1 N_scatD for gq 1 N_scat Assembling p q 4 2 Weight p q l alpha alphaxeta end end for p N_scatD 1 N_scat for gq 1 N_scat Assembling p q 5 3 Weight p q l alpha alphaxeta end end SCommon function A IntegralOperator 0O a M modes k Assembling Weight SSolving here direct density A B The post processing is then done by specifying to p diff how the density must be used for the first Mp obstacles a single layer potential is used while for the others a double layer potential is required The function can simply do that It just needs an array TypeOfOp of size N_scat x 2 such that 1 0 ifp lt Mp TypeOf0p p ie In a p diff script this means Apply the single layer potential multiplied by 1 for the first Mp part of the density and a double layer potential multiplied by 1 for the others sPreparing TypeOfOp TypeOfOp zeros N_scat 2 for p 1 N_scat if p lt N_scatD TypeOfOp p 1 else TypeOfOp p 2 end end sScattering angles theta_RCS 0 360 theta_RCS_rad theta_RCS 2 pi 360 SRadar Cross Section computation myRCS RCS O a M modes
10. either solve accurately the acoustic multiple scattering problem or to compute boundary integral operators for a collection of disks Indeed the diff toolbox can be used e as a solver of the multiple scattering problem e as a framework for more theoretical studies on boundary integral equations What does this user guide contain The first chapter recalls some well known properties of the boundary integral equations Even if this part is mostly theoretical some practical aspects of boundary integral equations are also provided for the impenetrable case and some examples of robust integral equations are given for the penetrable case Moreover some examples of usage are provided with the p diff toolbox Chapter B explains the theory on which the p diff toolbox is based More precisely the boundary integral operators are projected in the Fourier bases Their associated matrices are derived in addition to other physical quantities of interest such as the far or the near fields The p diff toolbox is detailed in chapter 3 each function is explained and examples are provided Simple example are provided in chapter 4 Let us note that the list of available functions is given in the appendix A in the alphabetic order and arranged by folder location in appendix B 8 Mathematical problem Reformulate the problem as an integral formulation Integral equation Pre processing geometry and right hand side b Assemble the mat
11. for p 1 M An example 2 6 is Np kap Cp where Cp weakly grows with kap A numerical study of the parameter Np is proposed in 2 6 where the following formula leads to a stable and accurate computation 1 Np kap 7 m 2V2rkape kap J 2 13 where is a small parameter related to the relative tolerance required in the iterative Krylov subspace solver used for solving the truncated linear system 2 14 see 2 6 The resulting linear system writes Lp U 2 14 where we introduced the block matrix L L i lt pq lt m the vectors p p icp lt m and U UP i lt p lt m as b1 b2 N LiM p U 21 b2 o L2M p U2 L oa pa WI U 2 15 LM LM2 p LMM p uM For p q 1 M the complex valued matrix L is of size 2N 1 x 2Nq 1 and its coefficients Lh are Lit LP for m Np Np n Ng Nq The complex valued components of the vector p ph n lt m lt n of size 2N 1 are the approximate Fourier coefficients p of p For the sake of clarity we keep on writing p p gt p for all m Np Np The complex valued vector UP U _n lt m lt n is composed of the 2N 1 Fourier coefficients of the trace of the incident wave on I i e U UP u p B2 LTY Ym Np Np If Niot TaN 1 denotes the total number of modes the size of the complex valued matrix L is then Mtot X Ntot More generally all th
12. ifq p The family Z m Z p 1 M also called Fourier or spectral basis is a Hilbert basis of L T for the usual scalar product pr 30 2 1 2 Integral operators integral equations for a cluster of circular cylinders In view of a numerical procedure u diff can use for example the weak formulation of the EFIE see Eq 1 29 page 22 in L T based on the Fourier basis Z Find p H 2 L such that for any p 1 M and m Z Lp PF ram ule Dh r20 Since u is assumed to be smooth enough typically and that T is then the scattered wavefield is also Q and the density p is at least in H 2 Therefore p can be developed in as p 3 gt ndh q 1 neZ and the weak form of the EFIE is Find the Fourier coefficients Pn C for q 1 M and n E Z such that Yp 1 M Vm EZ SE LEL 8 aq alr Bh am N This formulation can be written under the following matrix form Lp U where the infinite matrix representation L L i lt pq lt m and the infinite vectors p p icp lt m and U UP i lt p lt m are defined by blocks as Lil L2 LM p U 7 T21 22 fam 7 p 7 U L p U 2 1 pm TM o Mm aM iM with for any p q 1 M and m n Z Lea L R OF rary Ph pe and UP u lr Ph racry For the other integral formulations section or even for any other boundary condition the expressions of t
13. on a Matlab meshgrid and inside the obstacles e ExternalPotential Compute potentials single double or linear combination on a Matlab meshgrid and outside the obstacles PostProcessing NearField Functions e BlockPotential t Generic block potential matrix used to compute external or internal potentials e Options for potential computations are condensed here PostProcessing NearField Interface Internal potential of single layer potential only e ExternalSingleLayerPotential External potential of single layer potential only e ExternalDoubleLayerPotential External potential of double layer potential only e InternalDoubleLayerPotential Internal potential of double layer potential only PreProcessing Fourier Provide the number of mode to kept in the Fourier series PreProcessing Geometry Build a triangular lattice of disks Verify if the obstacles are satisfying the condition overlapping e RectangularLattice Build a rectangular lattice of disks Remove some disks e CreateRandomDisks Place randomly random disks in a box PreProcessing Incident Wave Full vector of a generic incident wave right hand side Block vector of a generic incident wave right hand side block 1 obstacle PreProcessing Incident Wave Block e Block vector 1 obstacle of a the right hand side of the trace of a point source wave Block vector 1 obstacle of a the r
14. potential only ut x Lt pt x Gre wpr y dy inQt u x 2 p x G x y y dy in r where te a 1 EA G FHS El yl corresponds to the Green function with respectively the wavenumber kt k and k kn note that k7 can differ from one obstacle to another Applying the transmission boundary conditions for the trace and normal derivative trace on I leads to the following couple system of boundary integral equations EN a F L pr 7 uihlp a 36 LENE AANA n j The upper index appearing on the boundary integral operators specifies which Green function is used in their expression For instance we have L p I G x y p y dy It can be proved that this integral equation is well posed if k is not an irregular frequency of the interior homogeneous Dirichlet problem The example BenchmarkPenetrable m located in Examples Benchmark solves the penetrable scattering problem by using this integral equation Actually the program has been written in the dielectric case not the acoustic one and thus 25 the transmission condition on the normal derivative is there slightly modified to allow a possible jump in the normal derivative Qnut spon ir Ayui p 0 where uo is the magnetic permeability of the medium and u py in Q is the magnetic permeability of the domain Q5 The associated integral equation is then obtained by modifying the quantity I See ies by in the matr
15. sin A x2 xs 2 x1 center 215 225 of a point source ais XS 1 and 25 XS 2 A point source wave is given by TAO klix xsl with x 1 2 and Xs 215 25 with HY the zeroth order Hankel function of the first kind k 1 x 1 wavenumber k in the vacuum k_int 1 x M wavenumbers in the obstacles k k_int p If k_int is a scalar then k5 k_int for all p 1 M M modes 1 x M vector of index of truncation of the Fourier series i e M_modes p Np Np 1 x 1 corresponds to the truncation index N in the Fourier series Nq 1 x 1 corresponds to the truncation index Nj in the Fourier series Integral operators see chapter 1 Index Letter diff abbreviation Operator 0 null operator 1 I Identity identity 2 L SingleLayer Lp fi G x y p y dy r 3 M DoubleLayer MA f Ony G x y A y dy r 4 N DnSingleLayer Np On J G x y p y dy r 5 D DnDoubleLayer DA a Ony G x y A y dy r 6 LIL PrecondDirichlet single scattering preconditioned trace of the single layer operator see 1 2 2 7 Dp PrecondNeumann single scattering preconditioned normal derivative trace of the double layer operator see 91 2 2 3 3 Pre processing The pre processing in diff consists in defining the right hand side or the incident wave and the geometry the obstacles The associated functions are located respectively in the folders mudiff PreProcessing Inci
16. tion for sound hard obstacles Here we only present this solution but the extension to other kinds of integral equations is direct The p diff script is close to the one developed for the Dirichlet problem only the two following functions must be modified is replaced by and the right hand side is now given by For the Neumann problem the preconditioned boundary integral equa tion is based on the double layer representation u Md D 1 DA D au oe Three unit disks 5 0 5 2 0 214 a y 1 1 sSet the parameters k 1 Swavenumber beta_inc 0 incident angle SFourier series truncation parameter M_modes FourierTruncation a k Min 1 Right hand side DnUincPrecond DnPlaneWavePrecond O a M_modes k beta_inc sMatrix of the system the two following lines are the same APrecond PrecondNeumann O a M_modes k SSolving here direct lambda APrecond DnUincPrecond The post processing is based on the double layer potential compared to Dirichlet the modifi cation is realized in the rcs function the last argument 1 0 is then replaced by 0 1 sScattering angles theta_RCS 0 360 theta_RCS_rad theta_RCS 2 pi 360 SRadar Cross Section associated with the double layer potential myRCS RCS O a M_modes k theta_RCS_rad lambda 0 1 plot theta_RCS myRCS k With the previous parameter we obtain the result shown on figure for an
17. 2Nq 1 1 Finally all these root vectors are stored in the three dimensional array Am 0 ifp q AR otherwise Am p q The global matrix A is thus stored through three different parts Az Am and Apr From the computer point of view as a three dimensional array has a fixed size in each direction the second and third dimensions of the arrays are both of size N_scat and the length in the first dimension is given by 53 3D array Length in the dimension 1 2 3 Ar max 2Np 1 N_scat N_scat Am ae ne 2Np 1 2N 1 1 N_scat N_scat qFP Ar max 2Np 1 N_scat N_scat With this constraint the vector Az p q can be larger than 2N 1 and thus A must be extracted A A 1 2N 1 p p The same occurs for Am and Ar Finally these three dimensional arrays are merged into a cell A representing the matrix A such that by using the Matlab notations a 1 Ar a 2 Am a 3 Ap Fast matrix vector products The matrix vector product between A and a vector X is divided into different elementary operations Let us consider Y AX and more particularly the pt component Yp X APX q By using the previous notations we have Yp ATX AT Air AR Xq g p Since AF At and Aji are diagonal matrices that are stored as vectors the corresponding matrix vector products are easy to compute The only difficulty concerns A4 Zq where Z4 ATX This c
18. BEMs For high frequency and or a large number of scatterers an adapted algorithm is moreover used that takes advantage of the the particular structure of the matrix block Toeplitz What is not p diff Even if u diff is designed to solve multiple scattering problem it is not a black box that only solves a particular problem i e p diff does not a priori fix the integral equation formulation to be used for solving a particular problem This theoretical aspect is rather let to the user However for users who are not fluent with boundary integral equations this user guide provides a rapid survey on boundary integral equations see chapter i Furthermore some examples of efficient solutions based on integral equations in classical cases Dirichlet Neumann boundary condition penetrable scatterers are detailled We expect that these elements will help a non specialist of boundary integral equation to use u diff for solving multiple scattering problems How to use u diff It is difficult to describe the way to use p diff without speaking about integral equations To simplify the presentation we provide a basic example Let us assume that a circular sound soft scatterer Q7 is placed within a homogeneous medium which is illuminated by an incident plane wave The mathematical problem is to compute the scattered field solution to the Helmholtz equation in the propagation domain A k u 0 in R Q u u on Q7 Sommerfeld s radi
19. FBWIE or if the far fields are not computed RCSSingleLayer RCS O a M modes k theta_RCS_rad rho 1 0 RCSBWIE RCS O a M modes k theta_RCS_rad psi_BW eta_BW 1 Near fields The scattered field can also be computed on a grid Even with a vectorization of the code this functionality remains slow and cpu consuming For the EFIE MFIE and CFIE only the single layer potential must be computed whereas for the Brakhage Werner integral equation the double layer must also be computed In all the case the computations are done in the exterior of the obstacles and ExternalPotent ial will do the job properly Building the grid XXmin 10 XXmax 10 YYmin 10 YYmax 10 le 0 1 XX XXmin lc XXmax YY YYmin lc YYmax 68 X Y meshgrid XxX YY Scattered field single layer potential ao 3 U ExternalPotential X Y O a M_modes k rho 1 0 OnBoundary 1 Scattered field Brakhage Werner integral equation UBWIE ExternalPotential X Y O a M modes k psi_BW eta_BW 1 OnBoundary 1 The total field wr u uinc can also be calculated thanks to which computes the incident wave on a grid even inside the obstacles on these point 0 value must be set using MaskMatrixObstacles 5 Incident wave UincOnMesh IncidentWaveOnGrid X Y k PlaneWave beta_inc sSet 0 to points in obs
20. Q5 p 1 M the single layer volume integral operator Y can be written as the sum of M operators 4 q 1 M defined by 4 HT Hig R Pq gt LaPa Vx E R Lapq x f G x y pq y dy q 21 Therefore the single layer potential can be decomposed as follows M Yp H V2 L Lp DEA with p4 Plr 1 27 q 1 By introducing the operators Lp q for p q 1 M defined on H 2 T by Voq HAR Lp aPa Lapa p that is Ypa H7 T Yx Ip Lpapa x G x y pq y dy 1 28 q then the EFIE can be written in the following matrix form L L 2 BEA iM p u r L21 L22 oe 2M p2 yine SE 1 29 gMt gp M 2 a LMM PM un n y The same procedure can be applied to the other volume operator and the other boundary integral operators M N and D to obtain ea HWE HP Lode f EEDA Mra MAT PPT MAg x f aG yA ely NPs OPT HPT Naglx f Angle youl APaly pea HAT HPT Da a f Onyx yA lAa q 1 2 2 Single scattering preconditioning Let us introduce the single scattering operator of the EFIE L as the diagonal part of the operator L defined by L OH 2 3 0 Ji 0 p 0 f 1 30 0 0 LMM Let k be a wavenumber k that is not an interior resonance k Fp Q7 Q implying that L and L are both invertible The single scattering preconditioner then simply consists in L L 1 1 0 ae 0 2 2 1 a 0 Her rte ose 0 0 0 ites a
21. incident plane wave of direction 0 and on figure 4 4 for a point source centered on 10 0 Only the case of the preconditioned integral equation has been shown 73 GMRES residual log Obstacles Radar Cross Section dB Radar cross section 4 8 a Obstacles History of convergence of the GMRES k 1 without restart O Precond 2 3 4 5 6 7 8 5 Iteration number c History of convergence of the GMRES Absolute value of the total field 1 8 16 14 12 1 o s os p 0 4 F 02 Jo s 6 4 2 o 2 4 6 8 io x e Absolute value of the total field 74 1 o 50 100 150 250 Angle of reception degree b Radar cross section Absolute value of the scattered field d Absolute value of the scattered field Real part of the total field 15 i l l 5 l 5 o 8 6 4 2 2 4 6 8 10 f Real part of the total field Figure 4 3 Neumann boundary value problem with an incident plane wave of direction 0 solved using only the single scattering preconditioned integral equation The three first figures show respectively the obstacles the radar cross section the history of convergence of the GMRES The three last d to f present the absolute value of the scattered field the absolute and real part of the total field Obstacles Radar cross section al all al a 1h 2 5 gt of 3 8 2 is al 3
22. of the eigenvalues of the Laplacian operator A for the homogeneous interior Dirichlet resp Neumann problem In the multiple scattering case that is when Q7 U OQ connected the following equalities hold true p is multiply M M Fp Q7 U Fo Q5 and Fy Q7 Fr 9 1 8 p 1 p 1 Throughout the paper Fpy Q7 denotes the set of all irregular frequencies Fon Q7 Fp Q7 J Fn 0 1 9 16 1 1 3 Direct integral equations Generalities This section details the way of deriving direct integral equations as described in 5 7 5 B3 This approach is nonstandard but has advantages that appear later in the user guide in section 1 2 21 The principle is to write the total field ur as a linear combination of a single and double layer potentials ur x LYp x W x u x Vx Ee Ot 1 10 where A p is now the unknown of the problem Thanks to proposition 1 2 such an expression ensures that both wr is solution of the Helmholtz equation in Qt and wp us is outgo ing Following 7 an integral equation is said to be direct when the densities A p have a physical meaning Indeed for these integral equations they are exactly the Cauchy data urlr Onur r However this is not a choice but a consequence of the construction of the integral equation In electromagnetic scattering direct and indirect integral equations are more often referred to as respectively field and source integral equations see e g Har
23. surface and far fields computations total and scattered exterior near field visualization 39 We now introduce the p diff toolbox based on Matlab by explaining the main predefined func tions and their relations with the previous mathematical derivations Section 3 2 gives a presen tation and introduces the notations Section 3 3 shows how to define the scattering configuration geometry and physical parameters Section presents the way the integral equations must be defined and solved Finally section 3 5 describes the data post processing The p diff toolbox is organized following the five subdirectories see also Figure 3 1 for a diagram of the directory arborescence e mudiff PreProcessing pre processing data functions incident wave and geometry section 3 3 e mudiff IntOperators functions for the four basic integral operators dense and sparse structure used in the definition of the integral equations to solve section 3 4 e mudiff PostProcessing post processing functions of the solution trace and normal derivative traces computation of the scattered total wavefield at some points of the spatial domain or on a grid far field and RCS section 3 5 e mudiff Common this directory includes functions that are used in p diff but which does not need to be known from the standard user point of view e mudiff Examples various scripts are presented for the user in standard configurations In addition
24. the p diff toolbox and first examples Contents oae be a kiuu a aia ad a aed WP ae Se ol se a 39 ar he Gee ae os xd a Soh eh a ad ca 40 3 3 PPe processing sss sso s amp we ee we ae Oe Se SS amp Se 42 3 3 1 Geometry creating the obstacles 42 3 3 2 Truncation of the Fourier series 2 0000 2 ee ee 46 3 3 3 Incident waves aaa aaa a 47 ee ee a a a he a ee a a a 49 BAL Generalities cis t oe 04 soa ee ea a ee ee we Rae a d 49 A ilable integr erator eri 50 did Aue MW ddd Altec ae ax ae add amp Oo Be eee wad 50 oo ts Seas eee aie feeds a eee AA eee A 52 3 5 Post Processing 0 0 2 eee ee ee eee ee eee ee 57 Show ade SS ee ak et a ee a A 62 3 1 Generalities An effort has been made to keep the same notations between the mathematical framework and the u diff toolbox Some explanations goal of the function and input output arguments about a function name_of_the_function can be obtained by typing the classical Matlab command window help name_of_the_function Because u diff includes all the integral operators that are needed in scattering traces and normal derivative traces of the single and double layer potentials a large class of scattering problems can be solved as long as the shape of the obstacles is circular Concerning the geo metrical configurations any deterministic or random distribution of disks is possible Finally p diff includes post processing facilities like e g
25. to the scattering problem A k u 0 in QF u e on Ip p 1 Mp nu nu onT p Mp 1 M u outgoing For this problem the preconditioned integral is not directly available We can apply a Combined Field Integral Equation for the mixed problem see equation 1 34 and written as I where the matrix A is given by 1 35 We have 1 a NP43 anl ifq lt Mp A p q p q a DP4 anM 4 if q gt Mp The matrix is particularly easy to build with p diff thanks to the frontal function IntegralOperator To this end two three dimensional arrays Assembling and Weight are built such that 4 2 ifq lt Mp and Weight p q 1 a an 5 3 if q4 gt Mp Sog Pp q n Assembling p q i The indices in Assembling corresponds to the indices of the boundary integral operators 2 LP4 3 M4 4 N 4 5 D41 The assembling process is realized by IntegralOperatoz Two Dirichlet obstacles unit diks OD 5 ape oy SIF aD 1 1 N_scatD length aD Two Neumann obstacles unit diks ON 5 5 5 5 aN 1 1 N_scatN length aN SA11 obstacles O OD ON a aD aN N_scat N_scatD N_scatN sSet the parameters k 1 Swavenumber beta_inc 0 incident angle SFourier series truncation parameter Dirichlet Neumann A11 M_modesD FourierTruncation aD k Min 1 M_modesN FourierTruncation aN k Min
26. 2 pi 360 R RCS O a M_modes k theta_RCS_rad rho_plus 1 0 The radar cross section is shown on figure 4 5 b Near field Let us now compute the near fields by first building the grid XX 20 0 1 20 YY 20 0 1 20 X Y meshgrid XX YY And second we compute the single layer potentials external and internal SSPotentials Ue ExternalPotential X Y O a M modes k rho_plus 1 0 Ui InternalPotential X Y O a M modes k_int rho_minus 1 0 OnBoundary 1 U Ue Ui The quantity U is the scattered field outside the obstacles and the total field inside in the obstacles The option OnBoundary is set to compute the potential also on the boundary in case there is some points of the grid on it To get the total field we need to add the incident wave to U outside the obstacles 79 UincOnMesh IncidentWaveOnGrid X Y k PointSource XS Matrix_Not_Obstacles MaskMatrixObstacles X Y O a 0 UincOnMesh UincOnMesh Matrix_Not_Obstacles U_tot U UincOnMesh The near fields are now fully computed and can be displayed Figure 4 5 c represents the absolute value of the total field and 4 5 d and 4 5 e its real and imaginary part 80 Obstacles Radar Cross Section T T o0000 FATF i OO0000 OO0000 WOOOO OO0O00 OO0O0O0 OO0OO0O O O O O O lt
27. 7 The unknown p is finally obtained through the solution of the direct integral equation Au 0 which can be written as AL p Au 1 14 Both the expression and the nature of the integral equation 1 14 depend on the boundary condition imposed to uy represented here by the operator A The next steps introduce three usual direct integral equations The proofs are not provided and can be found for example in 7 or 23 EFIE Electric Field Integral Equation For this first integral equation the operator A is the interior trace operator yg on I Thanks to the continuity on I of the single layer integral operator see equation 1 5 the boundary integral equation 1 14 becomes Lp u p 1 15 Due to theorem this first kind integral equation called Electric Field Integral Equation EFIE is well posed and equivalent to the scattering problem 1 2 except for Dirichlet irregular frequencies Proposition 1 5 If k g Fp Q then the function Zp u is solution to the scattering problem if and only if p is the solution of the EFIE 7 15 When k Fp Q7 the integral operator L is no longer bijective but is still one to one It can be shown that the kernel of the operator L is a subset of the kernel of the operator Z Consequently for every solution J of the EFIE the associated single layer potential Zp ute is still the solution of the scattering problem 1 2 MFIE Magnetic Field Integral Equation Ano
28. 9 otherwise 26 ap For the Neumann case the preconditioned integral equation EFIE reads as BD Baqui p The matrix representation D D of D D is then given by i IP if p q Vp q 1 M piz aq 2 7 dH 1 S 2 7 aJ otherwise ap As highlighted by Proposition there is no need to compute the preconditioned versions of the other integral equations as they lead to the same operator up to an invertible operators for BWIE To solve sound hard or sound soft scattering using the above integral equations appears to be a suitable choice 2 1 4 Projection of the incident waves in the Fourier basis To fully solve one of the integral equations EFIE MFIE CFIE or BWIE we need to compute the Fourier coefficients of the trace and normal derivative traces of the incident wave We give the results for both an incident plane wave and a pointwise source term Green s function For an incident plane wave the following proposition holds 2 Proposition 2 1 Let us assume that u cos 8 sin 8 and B 0 27 i e is an incident plane wave of direction B with B Yx R u x ekB x Then we have the following equalities UP up 8 JIm kap AU Onu p BF kin Jm kap p P L T L2 T with d forage Pare Let us consider now an incident wave emitted by a pointwise source located at s QF i e the wave u is the Green s function centered at s The Fourier coe
29. The p diff toolbox user guide Xavier ANTOINE and Bertrand THIERRY January 23 2015 version 0 9 Contents Copyright ntroduction to the user guide How to install y dif 1 Boundary integral equations a short survey 1 1 Standard integral equation formulations in acoustic scattering e tance RAR eas eee es os ee eae yo Sew A ehieas bok Gk GR Bow RBS Bw Ae ae Gee hae amp ner ee rr ances an ee Poe eee Gg eyes Re et GSEs Ge ee 1 1 6 Summary sius Hoe be OA Ay OR ee a ae HS DIESO 1 2 2 Single scattering preconditioning lt eoo e a a SOEREN The boundary value problem ooa a a OSOON pore e Fa a DO e ee ee E a ee ew RE ee Bee E 2 Multiple scattering by disks approximation method in u dif 2 1 Spectral formulation used in p difff 2 ee a 2 1 1 Notations and Fourier bases 0 0 000000 eee eee 2 1 2 Integral operators integral equations for a cluster of circular cylinders 2 1 3 Single scattering preconditioned integral equations 2 1 4 Projection of the incident waves in the Fourier basis 2 1 5 Near field evaluation 0 002 eee ee ee ee 2 1 6 Far field and Radar Cross Section RCS 220 0 000 1 4 1 2 2 Finite dimensional approximations and numerical solutions proposed in p diff 3 Description of the y diff toolbox and first examples 3 1 Generalities 3 2 Common argument and notations so a a e a 3 3 P
30. U has zero value outside the obstacles and by default on the boundaries too The argu ments are exactly the same as ExternalPotential Finally the common options for and are described below For the sake of clarity the examples are shown with only fExternalPotential even if they also apply for InternalPotential e ExternalPotential Verbosity VERBOSITY VERBOSITY is a scalar which when set to 0 stop displaying message default 1 e ExternalPotential OnBoundary ONBOUNDARY if ONBOUNDARY is set to 1 then the boundary is also covered beware of the jump rela tion default 0 59 Example 1 compute the single layer potential u p on a grid XX 10 0 1 10 yY 10 0 1 10 X Y meshgrid XxX YY U ExternalPotential X Y O a M modes k rho 1 0 Example 2 consider two obstacles and compute Lpi 0 5 x Lpz 1 5 x MA 1 5 X Mor2 inside the obstacle where the densities have already been computed XX 10 0 1 10 YY 1020 1 101 X Y meshgrid XX YY U InternalPotential X Y O a M modes k rho lambda 1 1 5 0 5 1 5 Remark 3 5 The near field functions ExternalPotential and InternalPotential are also interfaced by ExternalSingleLayerPotential InternalSingleLayerPotential for the single layer and by ExternalDoubleLayerPotential and InternalDoubleLayerPotential for the double layer potentia
31. Ur 0 0 A O 0 U gt 0 ge el eI 1 37 I _ Um 0 0 0 0o Ayu Now a second set of equations can be obtained thanks to the exterior part Following the same idea the trace and normal derivative trace of ut on I can be computed However the M contributions appear to be here l I M ra uplr ulr 5ut lr Do Mgg u le Lhal nur Vp 1 M i I Onutlr nu r Onu r Dt u r Nx Onu r 2 Pd Pd which can be rewritten as I Yp 1 M 505 agaj o with the following quantities o ull P M Lpg Up Onuplr ae es De an The operator Lyg is given by same for the three others LingPa bx yD ar ry q Finally this second set of equations can also be rewritten in the following matrix form I SAAT At AF Pa A 5 AT 1 2 13 1M uf ire 0 Ay 2A AR as us Uys 2 1 5 22 A23 2 M a nE DE i 1 38 Ut _ pine 0 M M Ant Am2 Am3 gt a Atm To obtain only one set of equations systems 1 37 and must be combined A first kind integral equation is obtained by doing 1 38 which is known to have a poor condition number but is robust and a second kind integral equation can be obtained by computing 1 37 1 38 which is well conditioned but can provide inaccurate solution These two integral equations are solved in the example function BenchmarkCalderon The matrix Aj is computed by the BlockCalderonProjector funct
32. a the right hand side of the trace of a point source wave PostProcessing NearField Functions Generic block potential matrix used to compute external or internal potentials BlockPrecondDirichlet IntOperators Dense Interface Block Block of the preconditioned integral operator sound soft case BLockP recondNeumann IntOperators Dense Interface Block Block of the preconditioned integral operator sound hard case BlockSingleLayer IntOperators Dense Interface Block Block of the single layer integral operator e BoundaryOfObstacles PostProcessing Geometry Extract the boundary of the obstacle CalderonProjector IntOperators Dense Interface Full Full dense matrix of the Calderon projector e CheckPlacement PreProcessing Geometry Verify if the obstacles are satisfying the condition overlapping e CreateRandomDisks PreProcessing Geometry Place randomly random disks in a box e dbesselh Common First derivative of first kind Hankel function e dbesselj Common First derivative of Bessel function e dbessely Common First derivative of Newton function DnDoubleLayer IntOperators Dense Interface Full Full dense matrix of the normal derivative of the double layer integral operator PreProcessing Incident Wave Full Full vector all obstacle of a the right hand side of the normal derivative of a plane wave e DnPlaneWavePrecond PreProcessing Incident Wav
33. acement function which can also be useful for a user who is placing the obstacles manually Removing disks The function aims to remove some disks of the geometrical configuration either disk by disk by row by column or by radius This can be useful for example to delete a row of disks Here is the syntax 45 O a RemoveDisk O_old a_old where O_old and a_old are the centers and radii of the current geometry Without optional argument the function has no effect The available arguments are e 0 a RemoveDisk X X1 X2 XN Remove all the disks with a center which has an abscissa equal to X1 X2 or XN e 0 a RemoveDisk Y Y1 Y2 YN Remove all the disks with a center which has an ordinate equal to Y1 Y2 or YN e 0 a RemoveDisk XY X1 Yl1 X2 Y2 XN YN Remove all the disks with a center with coordinates X1 Y1 X2 Y2 or XN YN e O a RemoveDisk Radius al a2 aN Remove all the disks with a radius equal to al a2 or aN e 0 a RemoveDisk Verbosity VERBOSITY set VERBOSITY to 0 to avoid display message to 1 to only show results and to gt 1 to see everything default e 0 a RemoveDisk Tol TOL Tolerance used for the conditional statement default 10719 Example 1 Remove all the obstacles that are either on the row of with abscissa equal to 1 or on the col
34. an however be achieved efficiently Indeed the discrete cross correlation product xcorr between A4 and Zq gives W xcorr Zg Adi where the bar denotes the complex conjugate of a complex number The result W is then extracted from W by Wa W 2Ng 1 2Ng 1 2Np The matrix vector product between A and X is then done in a fast way cross correlation is efficiently evaluated through the FFT Assembling the matrix The assembly process is very similar to the dense one the difference being the prefix Sp The block function is then called by SpBlockIntegralOperator Op ap Np Oq aq Nq Nmax k TypeOfOperator Weight The function returns three vectors corresponding to AV Af and AR The major difference here is that no linear combination of operators can be done during the assembly process for the sparse representation The quantity TypeOfOperator cannot hence be a vector but must be a scalar The assembly process of the global matrix can be done through 54 A SpIntegralOperator 0O a M modes k TypeOfOperator Weight The result is a Matlab ce11 of three components where each component is a three dimensional array The quantity TypeOfOperator can be a scalar or a matrix not a three dimensional array nor a vector If TypeOfOperator is a matrix then the block A 4 is assumed to be of type TypeOfOperator p q The linear combination of operators can still be done but must be specifi
35. and n i k Dense storage The operator I 1 a 5 anL is computed in p diff by using the IntegralOperator function the post processing remaining unchanged SCFIE alpha 0 5 eta i k sRight hand side Uinc PlaneWave 0O a M modes k beta_inc DnUinc DnPlaneWave O a M_modes k beta_inc BCFIE alpha xetax Uinc 1 alpha DnUinc 6 Assembling Matrix of the system the two following lines are the same ACFIE IntegralOperator O a M modes k 2 1 4 alphaxeta 0 5 l alpha l alpha Solving here direct rho ACFIE BCFIE Sparse storage The sparse storage version changes compared to the dense one the operators 2 N and L are computed separately and merged during the matrix vector products This is done in the SpMat Vec function SpL SpSingleLayer O a M_modes k SpN SpDnSingleLayer O a M_modes k SpA_MFIE SpAddIdentity SpN 0 5 M modes Solving and combining operators rho gmres X SpMatVec X M_modes SpL SpA_MFIE alphaxeta l alpha B_CFIE 4 1 5 The case of the single scattering preconditioned integral equation We strongly recommend to use the single scattering preconditioned version of the EFIE which is rigorously the same as the MFIE and CFIE and up to an invertible operator to any other boundary integral equation see Proposition 1 9 The EFIE version is available in diff and is represent
36. ation condition at infinity A possible integral representation of the scattered field u is the following u x G x y p y dy Yx ER 0Q where G is the Green function and p is the density solution of a boundary integral equation which can be e g the EFIE see chapter 1 Lo u with Lp x f G x y p y dy Vx ean dQ The problem is hence reduced to solving this boundary integral equation knowing p leads to the possible computation of u at any point x of the propagation domain Note that in addition the far field can be obtained easily thanks to its integral representation The theoretical part of the problem ends here and the computational part is handled now by p diff through the four following steps 1 Pre processing a Creating the obstacles b Building the right hand side u 2 Assembling the matrix of the integral operator s 3 Solving the resulting linear system 4 Post processing computing the far field near field at a point of the domain or on a grid or any other physical quantity of interest Most of the steps call some p diff functions except for the linear system solution which uses the already existing and optimized Matlab s functions for the direct backslash operator and iterative solvers GMRES This is summarized as a diagram in Figure i To whom diff is designed to The toolbox p diff is designed for any scientist who needs an easy to use and efficient way to
37. ce Far field of the double layer potential only e FarFieldSingleLayer PostProcessing FarField Interface Far field of the single layer potential only PreProcessing Fourier Provide the number of mode to kept in the Fourier series cet Potent ialOptions PostProcessing NearField Functions Options for potential computations are condensed here e HerglotzWave Examples TimeReversal FarField Common Compute an Herglotz wave linear combination of plane waves Ident it y IntOperators Dense Interface Full Full dense matrix of the identity operator PreProcessing Incident Wave Full vector of a generic incident wave right hand side e IncidentWaveOnGrid PostProcessing Incident Wave Compute incident wave on a Matlab meshgrid e IntegralOperator IntOperators Dense Generic integral operator dense matrix full e InternalDoubleLayerPotential PostProcessing NearField Interface Internal potential of double layer potential only InternalPotential PostProcessing Near Field Compute potentials single double or linear combination on a Matlab meshgrid and inside the obstacles 85 e InternalSingleLayerPotential PostProcessing NearField Interface Internal potential of single layer potential only e MaskMat rixObstacles PostProcessing Geometry Matrix with boolean values inside or outside obstacles P Lanewave PreProcessing Incident Wave Full Full vector al
38. ction IntegralOperator with two arguments the type of the operators for the identity operator and the double layer potential operator N see Table 3 2 and their associated weights 0 5 and 1 sRight hand side DnUinc DnPlaneWave O a M_modes k beta_inc Assembling sMatrix of the system the two following lines are the same A_MFIE IntegralOperator O a M modes k 1 4 0 5 1 Solving here direct rho A_MFIE DnUinc The post processing part is exactly the same as for the EFIE since the surface equation is based on the volume single layer integral representation Sparse storage The sparse storage version is almost the same as for the dense storage except for assembling the matrix and solving the linear system Indeed the matrices J and N cannot be com puted by the same function since the sparse function representations SpIntegralOperator and IntegralOperator cannot be summed together It is however possible to add the identity to a sparse operator thanks to SpAddIdentity SpN SpDnSingleLayer O a M_modes k S Add I 2 to N SpA_MFIE SpAddIdentity SpN 0 5 M modes rho gmres X SpMatVec X M_modes SpA_MFIE DnUinc 65 4 1 4 The case of the CFIE Let us now consider the well posed and well conditioned CFIE see also Eq 1 19 u Lp font 1 a G y p anu p 1 a dgu p Here we fix the parameters to a 0 5
39. ction 2 1 2 In addition low memory is only necessary when kap is large enough since the storage of the Toeplitz matrices can be optimized This resulting storage technique is called sparse representation in p diff in contrast with the dense full storage of the complex valued matrices Let us assume that ap a for 1 lt p lt M In terms of storage the dense version of a matrix requires to store about 4M ka coefficients assuming that N are fixed by formula 2 13 and r denotes the integer part of a real number r while the sparse storage needs about 4M ka complex valued coefficients In terms of computational time for solving the linear system the direct multi threaded gaussian solver included in Matlab leads to a cost that scales with O M ka For the preconditioned iterative Krylov subspace methods i e restarted GMRES the global cost is O M kalogs ka the converge rate depending on the physical situation and robustness of the preconditioner From these remarks we deduce that an iterative method can be an efficient alternative to a direct solver for large wavenumbers ka but also for large M We refer to 2 6 for a thorough computational study of the various numerical strategies A few examples in p diff are provided with the toolbox Finally the post processing formulas near and far fields quantities clearly inherits of the truncation procedure see sections and 2 1 6 37 38 Chapter 3 Description of
40. d Mathematics 204 2 440 451 2007 X Antoine C Chniti and K Ramdani On the numerical approximation of high frequency acoustic multiple scattering problems by circular cylinders J Comput Phys 227 3 1754 1771 2008 X Antoine and M Darbas Alternative integral equations for the iterative solution of acoustic scattering problems Quaterly J Mech Appl Math 1 58 107 128 2005 X Antoine and M Darbas Generalized combined field integral equations for the iterative solution of the three dimensional Helmholtz equation M2AN Math Model Numer Anal 1 41 147 167 2007 X Antoine and M Darbas Integral Equations and Iterative Schemes for Acoustic Scatter ing Problems to appear 2014 X Antoine K Ramdani and B Thierry Wide frequency band numerical approaches for multiple scattering problems by disks J Algorithms Comput Technol 6 2 241 259 2012 A Bendali and M Fares Boundary integral equations methods in acoustics in computa tional acoustic scattering In Computational Methods for Acoustics Problems pages 1 36 Saxe Coburg Publications 2008 T Betcke S N Chandler Wilde I G Graham S Langdon and M Lindner Condition number estimates for combined potential boundary integral operators in acoustics and their boundary element discretisation Numerical Methods for Partial Differential Equations 27 1 31 69 2011 S Borel R solution des quations int grales pour la diffraction d onde
41. d a point source centered on 10 0 the right hand side must be computed as XS 10 0 location of the source point source case sRight hand side 69 GMRES residual log Obstacles Radar Cross Section dB Radar cross section s EFIE e MFIE SpEFIE SpMFIE SpCFIE SpBW Precond Sparse Precond 4L fl 1 1 1 n i 1 f 1 1 f 6 4 2 o 2 4 6 1 0 250 300 350 x Angle of reception degree a Obstacles b Radar cross section History of convergence of the GMRES k 1 without restart Or 6 EFIE s MFIE Absolute value of the scattered field e CFE 10 1 e BW 2f Sparse EFIE P Sparse MFIE 09 Sparse CFIE Sparse BW 6 08 e Precond ar Sparse Precond i 0 7 2 06 6L gt 0 0 5 2 0 4 8L ad 0 3 6 0 2 104 8 0 1 10 o 12 f f 1 1 j 100 8a 6 4 2 o 2 4 6 8 10 0 5 10 15 20 25 x Iteration number c History of convergence of the GMRES Absolute value of the total field 1 6 14 1 2 1 0 8 0 8 0 4 0 2 o 8 6 4 2 o 2 4 6 8 10 s x e Absolute value of the total field 10 8 6 4 2 y o 2 4 6 8 10 A d Absolute value of the scattered field Real part of the total field 15 5 l 5 f o 8 6 4 2 2 4 6 8 10 f Real part of the total fie
42. dentWave and in mudiff PreProcessing Geometry 3 3 1 Geometry creating the obstacles The obstacles are stored in memory as a row vector a containing the radii of the disks and a 2 x N_scat matrix o with the centers of the disks 0 1 p a and 0 2 p y coordinates of the center of the pt disk These two arrays can be created either manually or by using some built in functions 42 Manually The disks can be created manually by simply creating the two variables and a containing respectively the coordinates of the disks and their radii For example for three obstacles placed at points 1 2 5 5 and 15 10 with respective radii 0 1 0 5 and 10 one gets O J l 2 2 7 L5 5y 10 a 0 1 0 5 10 Periodic placement Two built in functions are available with the toolbox to periodically place some disks with a rectangular or a triangular lattice as shown on Figure The two functions must be called as follows first for the rectangular lattice O RectangularLattice bx by Nx Ny OPTIONS and second for the triangular lattice O TriangularLattice bx by Nx Ny OPTIONS where the arguments are e bx distance separating two centers in the x direction be careful with the radii of the disks to avoid overlapping e by distance separating two row of obstacles in the x direction e Nx number of disks in a row e Ny number of rows For both functions the vector o
43. dition l a D an 51 m A a a u p am r l 1 24 This first kind integral equation is uniquely solvable for any wavenumber k and is equiv alent to the scattering problem BWIE the wave ur is represented as with S 7 4 0 Applying the boundary condition on I gives the BWIE 1 np I N p Beas 1 25 This first kind integral equation is also uniquely solvable for every k gt 0 and equivalent to the scattering problem 1 1 6 Summary The following table summarizes the main properties of the different integral equations for a Dirichlet respectively Neumann boundary condition Dirichlet resp Neumann boundary condition Integ Eq Nature Uniquely solvable for Correct physical solution EFIE 1 kind k Fp Q7 resp Fiy Q7 yes MFIE ont kind k g Fn Q7 resp Fip Q7 no CFIE 2 1 kind k gt 0 yes BWIE 2 1 kind k gt 0 yes 1 2 Multiple scattering case 1 2 1 A more explicit writing of the integral equation formulations The domain 27 is now supposed to be a collection of M disjoint bounded open sets Q of R p 1 M such that each domain R Qy is connected In the user guide single scattering designates scattering in a medium containing only one scatterer whereas multiple scattering is used for a medium containing more than one obstacle We focus on the multiple scattering case ie M gt 2 Since Q7 is composed of M disjoint obstacles
44. e Full Full vector all obstacle of a the right hand side of the normal derivative of a plane wave for the preconditioned problem of sound hard scattering PreProcessing Incident Wave Full Full vector all obstacle of a the right hand side of the normal derivative of a point source wave DnSingleLayer IntOperators Dense Interface Full Full dense matrix of the normal derivative of the single layer integral operator 84 e DORT_dielectric Examples TimeReversal FarField NonPenetrable DORT for penetrable obstacles far field e DORT_NonPenetrable Examples TimeReversal FarField NonPenetrable DORT for acoustic sound soft obstacles far field DoubleLayer IntOperators Dense Interface Full Full dense matrix of the double layer integral operator e ExternalDoubleLayerPotential PostProcessing NearField Interface External potential of double layer potential only e ExternalPotential PostProcessing Near Field Compute potentials single double or linear combination on a Matlab meshgrid and outside the obstacles e ExternalSingleLayerPotential PostProcessing NearField Interface External potential of single layer potential only e fangle Common Angle with horizontal axis PostProcessing FarField Generic far field computation from densities e FarField_to_RCs PostProcessing Far Field Radar Cross Section RCS from far field e FarFieldDoubleLayer PostProcessing FarField Interfa
45. e N24 k 2 4 LET An oe Car s ifp 4 imk ap i 3 dJ dH ifp q e Dea 2 5 ik2 eA w es E TVo aip ETa if p Ag where Spa is the transpose matrix of the separation matrix S 4 The integral equations involve the trace or normal derivative trace of the incident wavefield on T We have already introduced the infinite vector U of the coefficients of u p in the Fourier basis We then define similarly the infinite vector dU dU 1 lt p lt m Of the coefficients of the normal derivative trace Oyu p such that Vp 1 M meZ dU P aur OF ies Finally the density changes according to the integral equation and most particularly with respect to the boundary condition To keep the same notations as previously we introduce the densities and w used in the BWIE that are expanded in the Fourier basis as M M ee ke and Pa Ss X yD p l1meZ p l1meZ Finally we set A Neen and VW 1 lt pem where each block X Anez and UW WP mez is defined by Ym Z X AP and WP yP 32 2 1 3 Single scattering preconditioned integral equations The EFIE preconditioned by its single scattering component see Section 1 2 2 given by ELp Lou p can also be computed analytically in the Fourier bases Indeed let LIL be the matrix associ ated to the operator L L then Si P if p q 1 M L L PEET goin 2 6 a 4 fjr 1 Se T 7
46. e and bs is the normal derivative trace then B IncidentWave O a M modes k 3 4 1 2 In other words this builds the vector 3 4 can be translated as PointSource DnPointSource b lis Onui p inc a 1 uM HY klix xsl where with x 1 2 Remark 3 2 Remember that the resulting vector corresponds to the opposite of the trace or normal derivative trace 3 4 Integral operators 3 4 1 Generalities The functions defining the integral operators are available in the directory IntOperators which has the Dense and Sparse subdirectories for the dense matrix and sparse func tion representations of the four basic integral operators used in scattering i e L M N and given in their infinite dimensional operator versions by respectively 2 2 2 3 and 2 5 Preconditioned versions of the operators by their single scattering operators 24 are also defined Following Proposition 1 9 only L L EFIE Dirichlet and D D EFIE Neumann are provided since they are the only ones needed for the Dirichlet and Neumann problems Two different types of storage can be used within the p diff toolbox dense and sparse The assembly of the matrix is almost the same in both cases The way y diff is developed is close to the mathematics for the assembly process in the sense that the matrix A of an integral opera tor is built block by block A 4 For both storages two main functions
47. e boundary integral operators can be truncated according to this process Concerning the notations it is sufficient to formally omit the tilde symbol over the quantities involved in sections to 2 1 6 Since the four finite dimensional matrices L M N and D that respectively correspond to the four boundary integral operators L M N and D can be computed the linear systems that approximate the EFIE MFIE CFIE and BWIE can be stated For example the CFIE leads to with 0 lt a lt 1 and S n 0 lank 1 a G N p anU 1 a dU 2 16 Let us remark that the matrix obtained after discretization is always a linear combination of the four integral operators L M N D and the identity matrix I As a consequence for a given integral equation the resulting matrix is of size Niort X Ntot and has the same block structure as e g L see equation 2 15 The finite dimensional linear system or 2 16 is accurately solved in p diff by using the Matlab direct solver or a preconditioned Krylov subspace linear solver that uses fast matrix vector products based on Fast Fourier Transforms FFTs the choice of the linear algebra strategy direct vs iterative depending on the configuration with respect to kay and M The use of FFTs is made possible since the off diagonal blocks of the 36 integral operators can be written as the products of diagonal and Toeplitz matrices 2 6 see e g the matrices Spa in se
48. e located in the PostProcessing directory Far field The far field of the combination of the single and double layer potentials M 5 Np LpPp Yp Mprp p 1 is given by the equation 2 12 The corresponding u diff function is F FarField O a M modes k theta Density Weight where the parameters are e theta is the vector of receiving angles in radians e Density is either p A or both If Density is a column vector then Density 2 Density 1 and only one density is used both for the single and the double layer potentials contribution If Density has 2 columns then Density 1 is considered to be the single layer density and Density 2 the double layer density e Weight is the weight vector to apply to the M volume integral operators np and yp The quantity Weight is of size either 1 x 2 or 1x M The first column of Weight is applied to the single layer potentials and the second to the double layer potentials Weight p 1 p 57 and Weight p 2 p If Weight is a row of length 2 then Weight 1 n and Weight 2 p for all p the coefficient are the same for all obstacles With these notations the resulting far field is N_scat XO Weight p 1 Density 1 Weight p 2 Density 2 p 1 Example 1 the Dirichlet EFIE 1 15 the scattered field u reads as u p and the far field is then computed by F FarField 0O a M_modes k theta rho 1 0 E
49. e located in the directory called PostProcessing FarField Interface Near field inside and outside the obstacles The near field of a potential is given by equations and outside the obstacles and by 2 10 and inside In the same way p diff separates the outside and the inside computa tions Let us first explain how to compute the field outside the obstacles The user only needs to create a mesh by using e g meshgrid and launch ExternalPotential through the syntax U ExternalPotential X Y O a M_modes k Density Weight OPTIONS The resulting matrix U has zero value inside the obstacle By default the computation is not realized on the boundary of the obstacles this can be changed thanks to the options The arguments are e Density and Weight are exactly the same as for FarField e the quantities x and Y are the grid points which are created by the Matlab meshgrid function e the OPTIONS must be chosen from the options of the potential see below after the interior potential Let us consider now the computation inside the obstacle Q For the interior field the only contribution is assumed to be of the form L Pp Myp In other words the contributions of the other obstacles are not taken into account In the same way as for the external potential the computation is realized in p diff by the function call U InternalPotential X Y O a M_modes k Density Weight OPTIONS where
50. ed as u Lp E fT yin 66 Dense storage In p diff the quantity L hyine r is provided by P1anewaveP recond whereas L L is obtained with PrecondDirichlet The syntax for the dense version is then the following esad Right hand side UincPrecond PlaneWavePrecond O a M modes k beta_inc Matrix of the system the two following lines are the same APrecond PrecondDirichlet O a M modes k Solving here directly rho APrecond UincPrecond Pesa Sparse storage The sparse storage is here almost the same as for the dense version thanks to SpPrecondDirichlet SpPrecond SpPrecondDirichlet O a M modes k Solving and combining operators rho gmres X SpMatVec X M_modes SpPrecond UincPrecond 4 1 6 The case of the Brakhage Werner integral equation The Brakage Werner integral equation for the Dirichlet problem 1 21 reads as Dense storage As in the previous chapter we make here use of the common function IntegralOperator dared eta parameter eta_BW i k sRight hand side Uinc PlaneWave 0O a M modes k beta_inc sMatrix of the system the two following lines are the same A IntegralOperator 0O a M modes k 1 2 3 0 5 eta_BW 1 Solving here directly rho APrecond UincPrecond Dese Sparse storage For the sparse storage the three operators J L and N must be computed separately and mer
51. ed by Block SpBlockSingleLayer BlockSingleLayer The dense storage should be preferred for solving small scale problems most particularly low frequency problems when the memory storage and the CPU cost is not a problem where the sparse storage must be used when the limits are reached For the sparse storage which is detailed later the block matrices are stored through vectors and the matrix vector product is then fast when an iterative Krylov subspace solver is used This is essentially a suitable strategy when a sufficiently high number of modes N is required per obstacle i e for large enough wave numbers k The drawback is that some instabilities may arise in the numerical process if the truncations of the Fourier series are not done correctly The formula provides a stable result following 2 3 4 2 Available integral operators and numbering As for the incident wave for each operator there exists a function that builds the whole operator for the multiple scattering problem both for the dense and sparse storages The available operators are listed in Table 3 2 with their unique identifier integer their associated name in diff useful to get their interface functions and their definitions 3 4 3 Dense storage Building a block A An elementary block matrix in a global matrix can be created by using the following function BlockIintegralOperator Apq BlockIntegralOperator Op ap Np Oq aq Nq k Ty
52. ed in the matrix vector product see below Note that a function exists to add the weighted identity to a sparse operator see below Example 1 Creation of u diff sparse matrix of the single layer operator L L SpIntegralOperator O a M modes k 2 L is now a cell with three components Note the interface function can also be used m ll SpSingleLayer O a M_modes k Example 2 For two obstacles building the following matrix Min 0 5 x D192 1 5 x Doi 2 5 Xx N22 can be done with the commands TypeOfOp zeros 2 2 Weight zeros 2 2 TypeOfOp 1 1 3 Weight 1 1 1 block A_ 1 1 TypeOfOp 1 2 2 Weight 1 2 0 5 block A_ 1 2 TypeOfOp 2 1 5 Weight 2 1 1 5 tblock A_ 2 1 TypeOfOp 2 2 4 Weight 2 2 2 5 Sblock A_ 2 2 A SpIntegralOperator 0O a M modes k TypeOfOp Weight Assembling adding the identity It is possible to add the identity multiplied by a constant to a sparse operator A SpAddIdentity A alpha M_modes which simply returns A A al Example 1 If L is the sparse representation of the single layer operator then L SpIntegralOperator O a M modes k 2 alpha_L SpAddIdentity L 0 5 M modes computes 0 5 x L 55 Sparse matrix vector product A matrix vector product Y AX is done thanks to Y SpMatVec X M_modes ListOfOperators Weight
53. ere AP and ART are diagonal and called respectively the left and right parts and T Th with Te inem nq HO kbpq s is a Toeplitz matrix Tht Ta oa The idea is that the diagonal matrices can be stored as vectors containing the diagonal elements A matrix vector product between a diagonal matrix and a vector is then simply an element by element multiplication The Toeplitz matrix can also be stored in a compressed form as a vector and the matrix vector product is handled by a cross correlation based on the Fast Fourier Transform FFT The diagonals parts are stored in the three dimensional arrays Az and Ap L for Left part and R for Right part The left part Az also includes the diagonal part of the matrix diag A ifp q diag A otherwise 0 ifp 4q diag A s7 otherwise Ax 2 q l AR p q The matrix T of size 2N 1 x 2N 1 is then also compressed as a Toeplitz matrix following ty t2 t3 a t 2Nq 1 t 2Nq 141 ti t2 ok t 2Nq 1 1 TP t 2N 1 2 t 2Nq 1 1 ty ves t 2Nj 1 2 C 2Nq 1 2Nptl 1 2N4 1 2Np 1 2 F 2Ng 1 2Np 1 3 2NG41 2Np 1 1 Clearly the root vector containing the first row and the first column is enough to rebuild the matrix This root vector called A4 is of size 2N 1 2N 1 1 2N 2N 1 and such that t 2Nq 1 2Nptl1 1 t 2Nq 1 2Npt1 2 t 2N 1 2 API tN 1 1 ty tg t3 t
54. erface Block Block of the normal derivative of the double layer integral operator PreProcessing Incident Wave Block Block vector 1 obstacle of a the right hand side of the normal derivative of a plane wave PreProcessing Incident Wave Block Block vector 1 obstacle of a the right hand side of the normal derivative of a plane wave for the preconditioned problem of sound hard scattering 6 PreProcessing Incident Wave Block Block vector 1 obstacle of a the right hand side of the normal derivative of a point source wave IntOperators Dense Interface Block Block of the normal derivative of the single layer integral operator BlockDoubleLayer IntOperators Dense Interface Block Block of the double layer integral operator IntOperators Dense Interface Block Block of the identity integral operator e BlockIncidentWave PreProcessing Incident Wave Block vector of a generic incident wave right hand side block 1 obstacle 83 e BlockIntegralOperator IntOperators Dense Generic dense block of an integral operator PreProcessing Incident Wave Block Block vector 1 obstacle of a the right hand side of the trace of a plane wave PreProcessing Incident Wave Block Block vector 1 obstacle of a the right hand side of the trace of a plane wave for the preconditioned problem of sound soft scattering PreProcessing Incident Wave Block Block vector 1 obstacle of
55. es and My sound hard scatterers with Mp My M Without loss of generality let the sound soft obstacles Q be numbered for p 1 Mp and the sound hard ones for p Mp 1 Mp My The problem then reads as A k u 0 in QF u yine onT p 1 Mp nu O on Ip p Mp 1 Mp My wis outgoing to Q7 The field is then represented as a linear combination of single and double layer potentials M M um gt Ly pq D Maq q 1 q 1 where Z is given by 1 33 and by Mq HVT _ H R MO Mad WER Mg hq x A Ba O x y Ag y dy 83 A possible way to solve this problem is to apply the CFIE on both collections To derive the CFIE a fictitious field is introduced and the mixed condition 1 17 is applied on it Following the theory previously described this would lead to Ag 90 forqg 1 Mp Pq 9 for q Mp 1 Mp Mn 23 and the scattered field u is represented as Mp Mp Mn u gt gt Lipat XO Marq q 1 q Mp 1 To simplify let us assume now that there are only two obstacles the first one being sound soft and the second one being sound hard The extension to M various scatterers is straightfor ward The boundary condition applied to the interior fictitious field leads to the integral operator AfA AM Af Azma J where Ap represents the boundary condition 1 17 on I only Ap 1 g a Vp F ANYo p Applying the trace formulas leads to the system I 1 a 5 Nia anL 1
56. exist a block func tion BlockIntegralOperator and SpBlockIntegralOperator and a global matrix function Te eee This separation allows the user to either build a simple matrix for one operator or to con struct a more complex matrix where each block represents a different operator or a linear combination of them Let us also note that all the operators have interface functions located in Dense Interface or Sparse Interface folders For example SingleLayer is the interface function of to build the single layer operator and BlockSingleLayer the one for The same applies for the sparse storage with a prefix Sp 49 Int Op Identifier Function Definition for x T 0 zero operator null matrix I 1 Identity identity L 2 SingleLayer Lp x G x y p y dy Tr M 3 DoubleLayer MA x j Ony G x y A y dy r N 4 DnSingleLayer Np x On J G x y p y dy E D 5 DnDoubleLayer DX x a Ong GX y A V dy Tr Poy 6 PrecondDirichlet single layer preconditioned by its diagonal DD 7 PrecondNeumann double layer preconditioned by its diagonal Table 3 2 Available integral operators in y diff their unique identifier function name in terface and the mathematical definition The zero operator function does not have an inter face function and the sparse version is obtained by adding the prefix Sp SpDnDoubleLayer The block interface functions are also prefix
57. f is correctly installed on your computer and provides the expected results In the next chapter we present some standard examples of u diff scripts 62 Chapter 4 Simple examples of multiple scattering problems solved with u diff Contents seer oer ene sene 63 4 1 1 Pre processing 24444 wa See ae we a eA Wa ae bee eee 64 AAD Bat d Bb ded Ss ewe AM Ho ea a dd 64 Bias Be By epee A atenene 4 w 44 te 65 Ck ime he eae a ee ge 2 66 4 1 5 The case of the single scattering preconditioned integral equation 66 4 1 6 The case of the Brakhage Werner integral equation 67 4 1 7 Post processing gt s ee 68 ELS Results 4 ai a aaa a 9 209 8 E Bebe Boa we oe eA EE y R an 69 fa See ee hae Bea eree Me Ge goede 69 cos See OEN 73 Leveaewas 76 SiG eo Eo ee oe Re ae ER Oe we 78 4 4 1 Integral equation 2 2 20 e ee eee ee ee eee 78 4 4 2 A more complex geometry 0 0 a 78 4 4 3 Writing and solving the BIE using p difff 2 2 aa 78 4 4 4 Post processing 2 6 ee 79 The aim of this chapter is to provide some examples of multiple scattering problems solved by the p diff toolbox The impenetrable case with a Dirichlet a Neumann or a mixed of both boundary conditions set on the boundaries of the obstacles are fully treated For penetrable obstacles an example of the implementation of equation is provided 4 1 The Dirichlet boundary value problem Let us consider the scattering problem by a co
58. f radii must be built separately and manually For a set of unitary disks the following command is used a ones size O 2 The available options are same for RectangularLattice and TriangularLattice e RectangularLattice Origin Ostart places the first disk at Ostart posi tion Ostart being a 2x1 vector Default 0 0 e RectangularLattice Centered Ocenter center the collection on the point Ocenter Ocenter being a 2x1 vector Default none but replace Ostart value e RectangularLattice Direction DIR places the row in the increasing DIR 1 or decreasing x direction DIR 1 The default value is DIR 1 43 OOOO OO0O00 wn OOOO OWOVW00 1 OO000 10 0 0 0 0 a Rectangular lattice b Triangular lattice OO0O000 ak OO0O000 rae O00000 gk O00000 wb bl rb 000000 8 A Figure 3 2 Rectangular left and triangular right lattice with bx 3 by 4 Nx 5 Ny 6 Example Building a rectangular lattice of 3 x 4 unit disks Each disk is separated from the other by a distance equal to 1 5 so the centers are separated from a distance of 3 5 O RectangularLattice 3 5 3 5 3 4 a ones size 0 2 Random placement The toolbox p diff provides a function CreateRandomDisks to randomly place N_scat obstacles in a box xmin xmax x ymin ymax with a random radius In its simplest version the funct
59. f the identity operator Int Operators Sparse Interface Full e Sparse matrix of the single layer integral operator e SpDnDoubleLayer Sparse matrix of the normal derivative of the double layer integral operator e SpDnSingleLayer Sparse matrix of the normal derivative of the single layer integral operator e SpDoubleLayer Sparse matrix of the double layer integral operator Sparse matrix of the identity operator Sparse matrix of the preconditioned integral operator sound soft case Sparse matrix of the preconditioned integral operator sound hard case PostProcessing FarField Radar Cross Section RCS from far field Generic far field computation from densities Generic Radar Cross Section RCS computation from densities PostProcessing FarField Interface Radar Cross Section RCS of the double layer potential only Radar Cross Section RCS of the single layer potential only Far field of the single layer potential only e FarFieldDoubleLayer Far field of the double layer potential only PostProcessing Geometry e MaskMatrixObstacles Matrix with boolean values inside or outside obstacles e BoundaryOfObstacles Extract the boundary of the obstacle Display obstacles on figure PostProcessing Incident Wave 91 e IncidentWaveOnGrid Compute incident wave on a Matlab meshgrid PostProcessing NearField Compute potentials single double or linear combination
60. fficients of the trace and normal derivative trace of u on T are then given by the following proposition 23 33 nce is the Green s function Proposition 2 2 Lets Qt We assume that the incident wave u centered at s l TAO kllx sll The Fourier coefficients in of the trace and normal derivative trace of the incident wave on T are respectively given by Vx R s u x G x s 1 inc ima D UP u es Pe so T In kay HG kr p s Pin 8 and kagai inc ima T dU P Anu Irs BP saan hk P Jin Rap Hy hr p s Pn 5 2 1 5 Near field evaluation Outside the obstacles By using the Graf s addition theorem 18 23 we can compute the expression of the single and double layer potentials at a point x located in the propagation domain QT Proposition 2 3 Let p L T and p H 2 L be two densities admitting the following decompositions in the Fourier basis Z M M p X ph OP and SON A A p lmeZ p 1 mEZ Then for any point x in the domain of propagation Q the single layer potential reads SFe pP LBP x _ e ph i kap H krp x 8 x 2 8 p 1 mEZ p 1meZ and the double layer potential can be expressed as M M MN gt MBB 3 E MP ay HO kra Ra 29 p 1 mEZ p 1 mEZ Proposition 2 3 implies that for any x in QF ce PP Jn Kid X J kap EEC Ker xx BP x u x Lp x AX x a gt p 1meZ Inside the obstacles Sim
61. ficients of an incident wave either the trace of the normal derivative trace on one of the obstacles in the Fourier bases Its syntax is Bp BlockIncidentWave Op ap Np k TypeOfWave Param where TypeOfWave is a scalar value specifying the incident wave see table 3 1 and Param is the parameter of the wave angle of direction position of a point source The returned value Bp is a column vector of length 2Np 1 IncidentWave The function call is the following B IncidentWave 0O a M_modes k TypeOfWave Param The resulting vector B is of size yt 2M_modes p 1 The value Param is the same as for BLockIncidentWavel whereas TypeOfWave is a vector of size M where TypeOfWave p is the fixed choice for the block bp In other word the block bp is built by calling the function with the arguments BlockIncidentWave Op ap Np k TypeOfWave Param To simplify if TypeOfWave is a scalar value then it is considered as a vector with the same scalar value Example 1 Building a vector associated to the trace of an incident plane wave of direction beta_inc is done thanks to the following command the 1 argument refers to as PlaneWave 48 B IncidentWave O a M modes k 1 beta_inc or by using the interface function B PlaneWave O a M modes k beta_inc Example 2 For two obstacles and a point source centered at 1 2 if b is the trace of the wav
62. ged during each matrix vector products Note that in fact the operators J and N can be merged together as shown below using SpAddIdentity 67 SpL SpSingleLayer O a M modes k L SpM SpDoubleLayer O a M modes k M SpIminusM SpAddIdentity SpM 0 5 M modes 0 5I M Solving and combining operators psi_BW gmres X SpMatVec X M_modes SpL SpIminusM eta_BW 1 Uinc 4 1 7 Post processing The EFIE MFIE the CFIE and their preconditioned version share the same integral represen tation of the scattered field u as a single layer potential only u p Their post processing operations are the same For the Brakhage Werner integral equation the post processing is different since u M y Radar Cross Section RCS and far field Once the surface wavefield has been computed the far field can be calculated by the following diff commands here for a discretization of 0 27 with a step of 1 degree Post processing SScattering angles theta_RCS 0 360 theta_RCS_rad theta_RCS 2 pi 360 SFarfield for the single layer representation lt gt 1 0 FSingleLayer FarField O a M_modes k theta_RCS_rad rho 1 0 SFarfield for the BWIE lt gt eta_BW 1 FBWIE FarField 0O a M_modes k theta_RCS_rad psi_BW eta_BW 1 The RCS can also be computed simply by RCSSingleLayer FarFieldToRCS FSingleLayer RCSBWIE FarFieldToRCS
63. he only mathematical effort requires is that the user must of course choose the integral formulation that he wants to solve when the integral formulation is written then it can be solved by using p diff This chapter begins by presenting the potential theory and the four classical boundary integral operators with their main properties The case of the scattering by disks is then studied and the boundary integral operators are projected in the Fourier bases leading to infinite matrices but with some analytic expressions of their coefficients We then discuss the finite dimensional approximation The chapter concludes with the expressions of both the near and far fields and the projections of the right hand side onto the Fourier bases incident wave 13 1 1 Standard integral equation formulations in acoustic scatter ing We present a way to derive standard direct integral equations in the case of non penetrable obstacles Even if y diff can be used to solve the penetrable case studying the impenetrable case is a suitable way to introduce the boundary integral operators and their properties This section is strongly inspired by 5 7 15 23 1 1 1 Scattering problems Let us consider a homogeneous isotropic and non dissipative medium filling the whole space R and containing an open and bounded set Q7 possibly not connected but such that each component is itself simply connected Let I be its boundary and n the unit normal vector ou
64. he three boundary integral operators M N and D are needed Therefore to compute an integral equation we introduce the infinite matrices M MP lt pq lt M N NP 1 lt p q lt M and D lt p q lt M with the same block structure as L see equation 2 2 1 For p q 1 M the coefficients of the infinite matrices M EP NP41 and DP are derad for any indices m and n in Z by Mian Mh Pm rar Nmn NBh Pm rm and Dik DB Pm 121 For a numerical implementation we can explicitly compute 6 23 the matrix blocks LP 4 MP 4 N 4 and D involved in L M N and D for p q 1 M To this end we magine ihe infinite diagonal matrices JP dJ HP did dH with aaral terms for m Z Fram Im kap AD Pnm Jm kap Mem HG kap dl Fnm HG kap 31 In addition let P be the infinite identity matrix and for q p the infinite separation matrix S between the obstacles Q and 7 defined by SP4 meZ neZ and SP Drain Digg HE kbpg et ea Under these notations we rewrite the blocks LPa MP4 N and DP4 of the infinite matrices L M N and D under the matrix form for any p q 1 M irar Spip ifp q pets ER 2 2 Zn itp Aa lip _ kOe Fo aie lP eaj it p a Mer 2 2 2 3 ikt GpUgq z te Gray a if p q i ip 4 ikae e a if ya osi
65. ht 2 1 1 0 block A_ 2 1 TypeOfOp 2 2 1 4 Weight 2 2 0 5 1 Sblock A_ 2 2 A IntegralOperator 0O a M modes k TypeOfOp Weight Example 4 The Brakhage Werner Integral Equation BWIE for a Dirichlet problem 1 21 is solved by writing k 1 beta_inc pi eta i k Uinc PlaneWave O a M modes k beta_inc A IntegralOperator 0O a M modes k 1 2 3 0 5 eta_BW 1 psi A Uinc 3 4 4 Sparse storage Storing the matrices in a sparse way leads to a significant reduction in memory storage and the ability to get access to fast matrix vector product evaluations The linear system must then be solved by an iterative solver since the global matrix is no longer built Combining linearly the matrices is still possible but this is realized during the computation of the matrix vector products Indeed summing two matrices for two different integral operators is not guaranteed to keep the particular matrix structure If the problem involves at least two different integral operators they must hence be computed separately Compressed storage To explain how the matrix is stored by using the p diff sparse storage let us consider that A is a matrix corresponding to one of the four integral operators L M N or D The matrix A has the following special structure for p q 1 M and p q 52 e A is diagonal e A is full and can be written as AP ARTTP ARS wh
66. i 2max 1 R k leads to a reasonable condition number of the matrix of the linear system associated to the Brakhage Werner integral equation for a sufficiently high frequency Recent works have been developed on how to choose this parameter for much more general domains see for example 12 6 and 13 5 1 for the case of large k and 2 6 and 2 7 for the case of small frequency k Note also that according to Remark 2 24 these results apply to both Lgw and the CFIE operator since when a 1 2 these operators are adjoints up to a factor of 1 2 in the real L inner product Other generalizations of these equations when npw is an operator are available for example in 1 3 4 They provide a clear background concerning the generalization and improvement of both the CFIE and BWIE 19 1 1 5 Neumann boundary condition Consider now the scattering of a wave by a sound hard obstacle Neumann boundary condition Aut k u 0 in QF inc nu Onu on T u is outgoing The scattered field u is sought in the form of a linear combination between a single and a double layer potentials u x Lp x MANX xen To build the direct integral equations let us introduce the interior fictitious wave up 2 p MA wi defined in Q On T a boundary condition is hence enforced to the field Up such that it vanishes in Q In that case the jump relations 1 12 lead to p urlr p Onur p 0 The Neumann boundary condition the
67. ight hand side of the trace of a plane wave for the preconditioned problem of sound soft scattering 92 Block vector 1 obstacle of a the right hand side of the normal derivative of a plane wave Block vector 1 obstacle of a the right hand side of the trace of a plane wave Block vector 1 obstacle of a the right hand side of the normal derivative of a point source wave Block vector 1 obstacle of a the right hand side of the normal derivative of a plane wave for the preconditioned problem of sound hard scattering PreProcessing Incident Wave Full Full vector all obstacle of a the right hand side of the trace of a plane wave for the preconditioned problem of sound soft scattering Full vector all obstacle of a the right hand side of the trace of a plane wave Full vector all obstacle of a the right hand side of the normal derivative of a plane wave Full vector all obstacle of a the right hand side of the normal derivative of a point source wave Full vector all obstacle of a the right hand side of the trace of a point source Wave Full vector all obstacle of a the right hand side of the normal derivative of a plane wave for the preconditioned problem of sound hard scattering 93 94 Bibliography 1 F Alouges S Borel and D P Levadoux A stable well conditioned integral equation for electromagnetism scattering Journal of Computational and Applie
68. ilarly the potentials can be computed inside the obstacles which is useful for penetrable obstacles for instance In this case only the contribution of the current obstacle is taken into account u x LpPp Mp p Vx Qp Proposition 2 4 Let p L T and u H T be two densities admitting the following decompositions in the Fourier basis B M M p S_ S ph OP and A hoes p 1 mEZ p 1 mEZ 34 Then for any point x inside the obstacle Qp the single layer potential reads plo Tp LHR F D ph HE ky Jaler 10 meZ p l1meZ and the double layer potential can be expressed as 5 Y X MB x aS YN ITEC HO kap Jn hr p BF E 2 11 p l1meZ p 1meZ 2 1 6 Far field and Radar Cross Section RCS For computing the far field pattern let us recall that the scattered field u admits the following Helmholtz s integral representation u Z p MA where p and X are two unknown densities In the polar coordinates system r 0 and by using an asymptotic expansion of u when r 00 the following relation holds 14 etkr 1 V6 0 27 u r 7172 a z 0 aa 0 O Sa where ay and a y are the radiated far fields for the single and double layer potentials respec tively defined for any angle 8 of 0 27 by 1 iT ik av e yT 1 ik 8 n a f Q ik y X dr anl d r Iyl ye y y with 0 cos sin In addition the Radar Cross Section RCS is defined b
69. ing IntegralOperator creates a loop on all the obstacles p and q launches the following command Apq BlockIntegralOperator O p a p M modes p O q a q M modes q k Tpq Wpq and places Apq in the expected elementary block matrix of the global matrix The quantities Tpq and wWpq are given by TypeOfOperator and Weight such that weight is of the same size as TypeOfOperator and so Wpq follows the same rules as Tpq e If TypeOfOperator is a scalar then Tpg TypeOfOperator e If TypeOfOperator is a row or a column vector then Tpq is an array given by Tpq TypeOfOperator e If TypeOfOperator is a matrix then Tpq TypeOfOperator p q e If TypeOfOperator is a three dimensional array then Tpq is an array given by Tpq TypeOfOperator p q 51 Example 1 The single layer potential L is L IntegralOperator O a M_modes k 2 Example 2 The MFIE operator 0 5 x I N for a Dirichlet boundary condition 1 16 is A_MFIE IntegralOperator 0O a M modes k 1 4 0 5 1 Example 3 The following two obstacles operator is 0 5 x Thi Mii Li D21 0 5 x In N22 can be computed by using a three dimensional array and the null operator TypeOfOp zeros 2 2 2 Weight zeros 2 2 2 TypeOfOp 1 1 1 3 Weight 1 1 0 5 1 Sblock A_ 1 1 TypeOfOp 1 2 2 0 Weight 1 2 1 0 block A_ 1 2 TypeOfOp 2 1 5 0 Weig
70. ington and Mautz 16 9 or Borel 9 From now on the problem with unknown A p has only one equation given by the Dirichlet boundary condition on I To obtain a second equation a fictitious interior wave up defined in Q7 is introduced through u x Lo x M x u x Vx EeQ 1 11 Remark that on the one hand uzy is a solution of the Helmholtz equation in Q7 and on the other hand due to the trace relations 1 5 the couple of unknowns A p satisfies the well known jump relations A ur r urlr 1 12 p Onup r Onur r As the wave uy is fictitious it does not act on the solution ur of the scattering problem As a consequence the boundary condition on I imposed to up has no influence on up Let this constraint be represented by an operator A such that uy is the solution of the following interior problem ae in Q7 1 13 Aur 0 onr To build a direct integral equation the operator A is chosen such that the field up vanishes in Q Supposing that such an operator exists then the following equalities hold on the boundary T Fe 0 nuy r 0 Consequently and thanks to the Dirichlet boundary condition ur r 0 the jump relations 1 12 read as ea p Oyur r Therefore both the fictitious interior field up and the total field ur can be composed by a single layer potential only L p x ul x Vx 7 L p x ui x Vx EQ m Io B I g oe cad I 1
71. ion is called as O a CreateRandomDisks xmin xmax ymin ymax N_scat In that case createRandomDisks builds N_scat disks with unit radius in the box The function takes care to not overlap the disks Note that it is possible that the function does not succeed to place the obstacles e g if the user specifies too many obstacles in a box which is not large enough This is the reason why a security test has been set only 500 possible placements are allowed per disk The function comes along with a large set of optional arguments O a CreateRandomDisks xmin xmax ymin ymax N_scat amin amax dmim dmax O_avoid a_avoid dmin_avoid dmax_avoid where each additional argument is optional but the order must be kept amin must be set then amax etc and given by 44 Variable Type Default Description amin scalar 1 minimal random radius of the obstacles allowed amax scalar 1 maximal random radius of the obstacles allowed dmin scalar realmin minimal distance allowed between two obstacles not be tween the centers Setting lt 0 value will set dmin to realmin i e ignore it dmax scalar realmax maximal distance allowed between two obstacles not be tween the centers The maximal distance is efficiently reached Setting lt 0 value will set dmax to realmax i e ignore it o_avoia 2x N center of N hole s where the obs
72. ion and the whole matrix Aj 4 p q by CalderonProjector 27 28 Chapter 2 Multiple scattering by disks approximation method in p diff Contents ee e ae eee ee ere 29 2 1 1 Notations and Fourier bases oo oaa 2 020004 29 2 1 2 Integral operators integral equations for a cluster of circular cylinders 31 2 1 3 Single scattering preconditioned integral equations 33 2 1 4 Projection of the incident waves in the Fourier basis 33 BWR ay hy Bk Mech ar a a des ee Ra A Se aera 34 2 1 6 Far field and Radar Cross Section RCS o ooa 35 E a E hens 35 2 1 Spectral formulation used in u diff We consider that the scatterers are some circular cylinders In this case we can explicitly compute the boundary integral equations in some Fourier bases leading therefore to an efficient computational spectral method when used in conjunction with numerical linear algebra methods direct or iterative solvers 2 1 1 Notations and Fourier bases Let us consider an orthonormal system O Oz Ox We assume that the scattering obstacle Q is the union of M disks Q for p 1 M of radius ap and center Op We define T as the boundary of Q5 and by I Up 1 MIp the boundary of Q The unit normal vector n to Q is outgoing An illustration of the notations is reported in Figure For any p 1 M we introduce by as the vector between the center Op and the origin O b 00 bp bp Qp Angle Ozx1 bp
73. ions Sparse function add identity to a sparse operator e SpBlockDnDoubleLayer IntOperators Sparse Interface Block Sparse block of the normal derivative of the double layer integral operator e SpBlockDnSingleLayer IntOperators Sparse Interface Block Sparse block of the normal derivative of the single layer integral operator SpBlockDoubleLayer IntOperators Sparse Interface Block Sparse block of the double layer integral operator e SpBlockIdentity IntOperators Sparse Interface Block Sparse block of the identity operator 86 e SpBlockIntegralOperator IntOperators Sparse Generic sparse block of an integral operator SpBlockP recondDirichlet IntOperators Sparse Interface Block Sparse block of the preconditioned integral operator sound soft case 5SpBlockPrecondNeumanr IntOperators Sparse Interface Block Sparse block of the preconditioned integral operator sound hard case SpBlockS ingleLayer IntOperators Sparse Interface Block Sparse block of the single layer integral operator SpDnDoubleLayer IntOperators Sparse Interface Full Sparse matrix of the normal derivative of the double layer integral operator spDnsingleLayer IntOperators Sparse Interface Full Sparse matrix of the normal derivative of the single layer integral operator SpDoubleLayer IntOperators Sparse Interface Full Sparse matrix of the double layer integral operator SpIdent it
74. ix of 1 36 1 4 3 Calder n projectors Let us consider the Helmholtz representation of the interior field up u and the exterior one ul M utr 2 Anu p u umlo G wr Z nv r Vp 1 M where again the upper index specifies the interior or exterior Green function Note that the interior total field uy is equal to the interior field u However in the propagation domain the total field up is obtained by summing the scattered wave ut and the incident wave u The idea is to obtain two sets of M equations one taken from the interior and the other one from the exterior and next to combine them Let us first consider the M interior problems and take the trace and normal derivative trace of the interior fields see Eq for the trace formulae a I Nya z Up r 5 M up rp Lp nu Ir Vp 1 M I Baup lr Dy up rp 5 Np aup Ir Then by introducing the following quantities uz r Vp 1 M UL BoB P p ea the previous set of equations can be rewritten as I Yp 1 M 452 5 0 where Azs My Ly Dpp Npp The operator Lpp are the operators L with the interior Green functions same for the three others LppPp G klix yll Pp arp Iry P 26 Note that the Calder n projector P defined by I Pp E 2 F Ap is here hidden This set of equations can be rewritten in the following matrix form I 0 O ie 0 5 1 1
75. k theta_RCS_rad density TypeOfOp plot theta_RCS myRCS k 1 m TT 4 4 Penetrable case 4 4 1 Integral equation Let us now consider the case of penetrable obstacles with a frequency k possibly different in each scatterer Q for p 1 M The problem then reads as A kt ut 0 in Qt A k7 uT 0 in Q7 Cy ue onT nut nu dyguvi onT ut outgoing where k o kp Solving this problem can be done thanks to the integral equation presented in section where the interior total field u equal also to the total field and the exterior scattered field ut are represented as a single layer potential ut x Lt p x Gt e y or y dy in Ot u x 7 p x aa y e y dy inQ The associated integral equation 1 36 is recalled to be _ rs i L pt 7 u a a ar 2 4 4 2 A more complex geometry As already explained in 3 3 1 y diff also proposes some tools to place the obstacles and in a particular order leading to original geometry as the one shown previously on figure and treated here see 4 5 a The geometry is composed by a periodic placement of 11 x 11 unit disks separated by a distance equal to 1 with an empty middle line and middle column leading to in fact 10 x 10 100 obstacles Below is recalled the code to obtain this geometry bx 3 by 3 Nx 11 Ny 11 O RectangularLattice bx by Nx Ny Centered 0
76. l These four ready to use functions are located in the following folder PostProcessing NearField Interface Incident wave The incident waves can also be computed on a grid by using Uinc IncidentWaveOnGrid X Y k TypeOfWave Param where like for the Incident Wave function from the pre processing part see 3 3 3 we have e X and Y are two matrices coming from the meshgrid Matlab function e Param is either the incidence angle of a scalar plane wave or the location of a point source 2 x 1 vector e TypeOfWave is either a char or a scalar value Example computation of a plane wave with an incidence angle equal to beta_inc 7 beta_inc pi XX 10 0 1 10 YY 1020 1 10 X Y meshgrid XxX YY Uinc IncidentWaveOnGrid X Y k PlaneWave beta_inc The last line could also be switch PlaneWave with its identifier 1 Uinc IncidentWaveOnGrid X Y k 1 beta_inc 60 Geometry drawing the obstacles and creating mask matrices The p diff toolbox also proposes some functions for the post processing of the obstacles Drawing the circular cylinders can be done by calling the following function PlotCircles O a fig_index OPTIONS where fig_index is the Figure handle and the OPTIONS are e PlotCirclesOnFigure Color COLOR apply the color COLOR to lines same as the plot function e PlotCirclesOnFigure
77. l obstacle of a the right hand side of the trace of a plane wave PreProcessing Incident Wave Full Full vector all obstacle of a the right hand side of the trace of a plane wave for the preconditioned problem of sound soft scattering PostProcessing Geometry Display obstacles on figure PreProcessing Incident Wave Full Full vector all obstacle of a the right hand side of the trace of a point source wave PrecondDirichlet IntOperators Dense Interface Full Full dense matrix of the preconditioned integral operator sound soft case PrecondNeumann IntOperators Dense Interface Full Full dense matrix of the preconditioned integral operator sound hard case RCs PostProcessing FarField Generic Radar Cross Section RCS computation from densities PostProcessing FarField Interface Radar Cross Section RCS of the double layer potential only PostProcessing FarField Interface Radar Cross Section RCS of the single layer potential only e RectangularLattice PreProcessing Geometry Build a rectangular lattice of disks PreProcessing Geometry Remove some disks ee e repeat_horiz Common Copy paste a row vector to build a matrix e repeat_vert Common Copy paste a column vector to build a matrix SingleLayer IntOperators Dense Interface Full Full dense matrix of the single layer integral operator spaddident ity IntOperators Sparse Funct
78. ld 10 8 6 4 2 y 0 2 4 6 8 10 A Figure 4 1 Different results for a Dirichlet problem and an incident plane wave solved using p diff The first figure shows the three obstacles and the second one the radar cross section obtained for the different integral equation all superimposed The next figure c represents the GMRES history of convergence of the different integral equations The figure d shows the absolute value of the scattered field and the last two e and f represent respectively the absolute value and the real part of the total field 70 Uinc PointSource O a M_modes k XS The other change appears in the post processing when computing the total field on a grid Indeed the incident wave is different and IncidentWaveOnGrid must get the right argument UincOnMesh IncidentWaveOnGrid X Y k PointSource XS The rest of the code remains unchanged and the results obtained with a point source are shown on figure 71 Obstacles Radar cross section al 3b g SpEFIE S SpMFIE SpCFIE al p gt a 1H Zz E gt gt of o 8 k 6 T at 5 G ba a Oo z B al al 6 4 2 0 2 4 6 150 200 250 300 350 x Angle of reception degree a Obstacles b Radar cross section History of convergence of the GMRES k 1 without restart e EFIE e MFIE Absolute
79. ld and combined field solution for conducting bodies of revolution Archiv Elektronik und Uebertragungstechnik 4 32 157 164 1978 R Kress and W T Spassov On the condition number of boundary integral operators for the exterior Dirichlet problem for the Helmholtz equation Numer Math 42 1 77 95 1983 P A Martin Multiple Scattering Interaction of Time Harmonic Waves with N Obstacles volume 107 of Encyclopedia of Mathematics and its Applications Cambridge University Press Cambridge 2006 J Mautz and R Harrington A combined source solution for radiation and scattering from a perfectly conducting body Antennas and Propagation IEEE Transactions on 27 4 445 454 jul 1979 W C H McLean Strongly elliptic systems and boundary integral equations Cambridge University Press 2000 J C N d lec Acoustic and Electromagnetic Equations Integral Representations for Har monic Problems volume 144 of Applied Mathematical Sciences Springer Verlag New York 2001 Y Saad and M Schultz GMRES a generalized minimal residual algorithm for solving nonsymmetric linear systems SIAM J Sci Statist Comput 7 3 856 869 1986 B Thierry Analyse et Simulations Num riques du Retournement Temporel et de la Diffrac tion Multiple Nancy University Th se de Doctorat 2011 B Thierry A remark on the single scattering preconditioner applied to boundary integral equations Journal of Mathematical Analysis and Applicati
80. ld solve the multiple scattering problem based on classical boundary integral equations as described in chapter I and displays the Radar Cross Section and the history of GMRES for the various boundary integral equations as Figure 4 Other tests can be launched to check that all evaluations of the integral operators works correctly BenchmarkNeumann BenchmarkPenetrable BenchmarkCalderon 5 If everything went right congratulation Your p diff installation is working 11 Obstacles a Obstacles Radar cross section Radar Cross Section dB 1 e EFIE o MFIE o CFIE e BW SpCFIE SpBW s Precond Sparse Precond 2 Angle of reception degree 250 c Radar cross section Figure 2 4 of the different figures that pop up after launching BenchmarkDirichlet 24 4b GMRES residual log fl b 84 Absolute value of the total field b Absolute value of the total field History of convergence of the GMRES k 1 without restart Sparse EFIE Sparse MFIE Sparse CFIE Sparse BW Precond Sparse Precond 2 4 6 8 10 12 14 16 18 20 Iteration number d History of convergence of the GMRES The first figure shows the two obstacles and the second one the total field The two others figures c and d represent respectively
81. llection of sound soft obstacles A k u 0 in QF u u on u outgoing 63 with Q7 U Q We propose to solve this problem through various integral equations the EFIE 13 the MFIE 1 16 the CFIE 1 19 and the single scattering preconditioned integral equations 1 32 We show how to use both the full and sparse storages of the matrices Before starting we recommend to use the single scattering preconditioned integral equation as presented in Indeed the resulting system is well posed and is well conditioned leading to an efficient solution by a Krylov subspace iterative solver 4 1 1 Pre processing Let us first consider a collection of three sound soft unit circular cylinders The wavenumber is k 27 and the direction of incidence of the wave is 8 0 degree The resulting p diff pre processing code for setting these parameters is then 6 Pre processing Three unit disks oe O 5 0 5 2 0 2 a 1 1 1 sSet the parameters k 1 Swavenumber beta_inc 0 incident angl plane wave case SFourier series truncation parameter M_modes FourierTruncation a k Min 1 For each integral equation we now present the assembly process the computation of the solution and finally the post processing of the computed wave fields The common pre processing part is the one described above All the functionalities presented here are also available in the file BenchmarkDirichlet m which is l
82. n makes the density p vanishing The interior up and exterior ur fields are obtained through a double layer potential uz x MA x u x EQ ur M x u x Vx Qt As previously the boundary condition applied to uy is the integral equation to solve Here are listed the four integral equations obtained in addition with their properties e EFIE applying a homogeneous Neumann boundary condition to the wave up Onur r 0 1 22 leads to the EFIE for the Neumann boundary condition DA Ayu p This first kind integral equation is well posed and equivalent to the scattering problem as long as k is not an irregular frequency for the interior homogeneous Neumann problem If k Fy Q7 however then after reconstruction the solution obtained by the EFIE is correct e MFIE applying a Dirichlet boundary condition to the field uz uz r 0 leads to the MFIE for the Neumann boundary condition 1 p 3 M A u p 1 23 This second kind Fredholm integral equation the operator M is compact is well posed It is equivalent to the scattering problem if k is not an interior resonance for the Dirichlet boundary value problem In that case the solution obtained from the MFIE is not correct and must not be used for a practical computation 20 CFIE applying the mixed boundary condition 1 17 to the wave up 1 a nur r anuyr r 0 with O0 lt a lt l et S n 4 0 leads to the CFIE for a Neumann boundary con
83. nUinc DnPlaneWave O a M_modes k beta_inc rho gmres X SpMatVec X M_modes SpI SpN 0 5 1 DnUinc or by using the SpAddIdentity function 56 SpAMFIE SpDnSingleLayer O a M modes k SpAMFIE SpAddIdentity SpAMFIE 0 5 M modes DnUinc DnPlaneWave O a M_modes k beta_inc rho gmres X SpMatVec X M_modes SpAMFIE DnUinc Example 3 for the Dirichlet BWIE 1 21 The operators L and M must be computed separately and 0 57 can be added to M directly k 1 beta_inc pi eta i k Uinc PlaneWave O a M modes k beta_inc SpL SpSingleLayer O a M_modes k SpM SpDoubleLayer O a M_modes k SpIminusM SpAddIdentity SpM 0 5 psi gmres X SpMatVec X M_modes SpL SpIminusM eta_BW 1 Uinc 3 5 Post Processing Now that the system has been solved the next step is to display some results The u diff toolbox proposes some post processing features such as computing the far field fast or the near field of the wave on a grid Matlab meshgrid or only at some points In addition other possibilities are offered such as drawing the disks on a figure or the incident field Both the near and far field computations can be done for a linear combination of a single and double layer potentials Lp MX x only one of them e g Zp or with the same density Z W w with a single command line All post processing functions ar
84. ocated in the Examples Benchmark folder 4 1 2 The case of the EFIE This integral formulation reads as u Lp Lp u r where the first line is the integral equation representation of the exterior wavefield u and the second one is the surface integral equation to solve Dense storage In p diff the surface single layer operator L and the incident plane wave field u p are prede fined quantities If the full storage of the integral equation is used then the direct solution of the resulting linear system can be obtained by the standard backslash Matlab operator sRight hand side plane wave Uinc PlaneWave 0O a M modes k beta_inc Assembling sMatrix of the system the two following lines are the same L SingleLayer O a M modes k Solving here direct rho L Uinc oe 64 Sparse storage For the sparse storage version only the assembly process of the single layer matrix and the system solution need to be modified as follows sMatrix of the system the two following lines are the same SpL SpSingleLayer 0O a M_modes k Solving here direct rho gmres X SpMatVec X M_modes SpL Uinc 4 1 3 The case of the MFIE The resolution of the scattering problem by the MFIE 1 16 leads to the integral equation representations u Lp I Lan anp Dense storage The MFIE operator I N GA can be computed thanks to the frontal fun
85. of the normal derivative of the double layer integral operator Full dense matrix of the normal derivative of the single layer integral operator Full dense matrix of the double layer integral operator Full dense matrix of the single layer integral operator Full dense matrix of the Calderon projector Full dense matrix of the preconditioned integral operator sound hard case Full dense matrix of the preconditioned integral operator sound soft case e Identity Full dense matrix of the identity operator IntOperators Sparse e SpBlockIntegralOperator Generic sparse block of an integral operator e SpIntegralOperator Generic integral operator sparse matrix full IntOperators Sparse Functions e Sparse function add identity to a sparse operator Sparse function sparse matrix vector product possibly multiples e Sparse function sparse matrix only one vector product Int Operators Sparse Interface Block 90 Sparse block of the double layer integral operator Sparse block of the preconditioned integral operator sound hard case Sparse block of the preconditioned integral operator sound soft case e SpBlockDnSingleLayer Sparse block of the normal derivative of the single layer inte gral operator e SpBlockDnDoubleLayer Sparse block of the normal derivative of the double layer integral operator Sparse block of the single layer integral operator Sparse block o
86. ometry 00 00 0 ee ee 4 4 3 Writing and solving the BIE using p difff 2 4 4 4 Post processing 2 ee A List of the diff functions alphabetical order B List of the diff functions ordering by folder name Copyright Copyright 2014 Xavier Antoine Bertrand Thierry Universit de Lorraine Institut Elie Cartan de Lorraine UMR CNRS 7502 F 54506 Vandoeuvre l s Nancy Cedex FRANCE Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies Introduction to the user guide What is p diff The toolbox p diff has been developed for solving two dimensional acoustic multiple scatter ing problems by disks based on boundary integral equations The toolbox is a set of Matlab functions that compute accurately and efficiently the standard integral operators with in ad dition pre and post processing facilities When the boundary integral formulation has been written mathematically the coding part can be easily done in p diff Solving the linear system is realized within the Matlab framework No particular computational skill is needed except basic notions in Matlab The boundary integral operators are spectrally discretized in Fourier basis thanks to the circular shape of the obstacles This leads to an accurate solution and a smaller size linear systems to solve compared to usual Boundary Element Methods
87. ons 413 1 212 228 2014 96
88. peOfOperator Weight The Weight argument is optional and set to 1 by default The quantity TypeOfOperator speci fies the integral operator to compute thanks to the numbering of Table 3 2 If TypeOfOperator 50 is a scalar e g 2 then the resulting matrix Apq is the elementary matrix of the associated operator e g L 41 If TtypeofOperator is a row e g 1 3 the sum of the two operators is computed e g 1 4 M Finally the Weight quantity of the same size as TypeOfOperator is the constant that is used to multiply the block and hence N API 5 Weight f Operator 1 where Operator f is one of the integral operators Example 1 Build the single layer block L Apq BlockIntegralOperator Op ap Np Oq aq Nq k 1 Example 2 Build the whole matrix 0 5 x I 4 NP14 appearing in the MFIE 1 16 Apq BlockIntegralOperator Op ap Np Oq aq Nq k 1 4 0 5 1 Example 3 Compute the sum of the four blocks 0 5 x L 4 1 5 x MP4 2 5 x N 443 5 x DP Apgq BlockIntegralOperator Op ap Np Oq aq Nq k 2 3 4 5 Or Sr L 5 2 5 345 7 Assembling the matrix A Now that the construction of an elementary block matrix is well understood building a global matrix is easy thanks to the common function A IntegralOperator 0O a M modes k TypeOfOperator Weight As for the quantity Weight is optional and set to 1 by default Roughly speak
89. penetrable scattering using Calderon pro jectors e BenchmarkDirichlet Example of solution of sound soft scattering using different in tegral equations e BenchmarkNeumann Example of solution of sound hard scattering using different inte gral equations e BenchmarkPenetrable Example of solution of penetrable scattering using single layer potential Examples TimeReversal FarField Common e HerglotzWave Compute an Herglotz wave linear combination of plane waves e TimeReversalOperator Time reversal matrix in acoustic and far field context Examples TimeReversal FarField NonPenetrable e DORT_dielectric DORT for penetrable obstacles far field e DORT_NonPenetrable DORT for acoustic sound soft obstacles far field Int Operators Dense 89 e BlockIntegralOperator Generic dense block of an integral operator e IntegralOperator Generic integral operator dense matrix full Int Operators Dense Interface Block e Block of the normal derivative of the single layer integral oper ator Block of the double layer integral operator Block of the identity integral operator Block of the preconditioned integral operator sound soft case Block of the preconditioned integral operator sound hard case Block of the single layer integral operator Block of the normal derivative of the double layer integral oper ator Calderon projector block IntOperators Dense Interface Full e Full dense matrix
90. r by M HYT HERAT A MA Yx R T W x f On G x y A y dP y T where the spaces H 12 T H R T Hi R T are the usual Sobolev spaces and the two dimensional Green function G is given by i THP klix yl 1 3 Vx y ER x y G x y The function HP is the first kind Hankel function of order zero The way the integral operators on T must be understood is through the duality product between the Sobolev space H 2 T and its dual space H 2 P However when the data u and T are smooth enough then the scattered field u is also regular and the dual product can be identified to the non hermitian inner product on L I Fo me FI a x This identification is considered throughout this paper Let us now define the trace y and the normal derivative trace yf operators see 13 Appendix AJ where the plus minus signs specify whether the trace is taken from the inside of Q Q7 First the trace operators yf H1 Q H L are defined so that if v C QE then xo li a ae ue for almost every x I By introducing the space H Q F A v H QF Av L 0 and the linear operators y H 2 gt H1 Q such that 7729 4 for all y H T the normal traces yF H Q A gt H L can be defined 13 Equation A 28 Vu H1 0 A Vy HY PL tv p H 1 2 1 H1 2 r 7 i i Av x w x d
91. r any frequency k Proposition 1 7 For any k gt 0 and for any complex valued numbers a and 7 satisfying condition the function Zp u is the solution of the scattering problem if and only if p ts the solution of the CFIE 7 2 1 1 4 Brakhage Werner indirect integral equation The indirect Brakhage Werner Integral Equation BWIE is derived now Let us first start by introducing some notations The total field wr is here sought as a linear combination of a single and a double layer potentials applied to the density y H L up u Lowi where the operator gw with parameter npw is given by Vx R T Yaw x mw Myx 1 20 f 20 G0x mwa y Wy ro 1 20 with S n 0 The integral equation is obtained by applying the exterior trace y to ur Indeed the Dirichlet boundary condition ant ur 0 and the trace relations 1 5 directly give the Brakhage Werner integral equation with 7 as unknown Law u p 1 21 with i Lgw nz M 51 This second kind integral equation does not suffer from irregular frequency 0 Proposition 1 8 For all k gt 0 the quantity Zgyw u is the solution to if and only if w is the solution of the Brakhage Werner integral equation 7 21 A numerical study concerning the optimal choice of the parameter npw see relation 1 20 is proposed in I7 in the case of a single spherical or circular obstacle of radius R For a Dirichlet boundary condition the choice ngw
92. re processing 3 3 1 Geometry creating the obstacles 2 02 0004 3 3 2 Truncation of the Fourier series 0 0 0000 ee eee 3 3 3 Incident waves 11 13 14 14 14 17 19 20 21 21 21 22 23 24 24 25 26 31 3 4 Integral operators ooo ee 3 4 1 Generalities 2 ee 3 4 2 Available integral operators and numbering 3 4 3 Densestoragel sis eee eR Re DAS Gs Bae we Se ER Goh eo 3 44 Sparse StOragel yos see Gs sek ee eek ee eo Se Boe e ee Get ES 350 Post Processitigj sog sogo foe ec Gok Goa RO ee ee Oe Ree Eee dee 3 6 Examples available in u difff 2 ee ee 4 Simple examples of multiple scattering problems solved with p difi 4 1 The Dirichlet boundary value problem 202 0004 4 1 1 Pre processing m pi ets s a ee Re Soe ee Re Re ee Ea ke eG sth Wie Te ee ee eo ee Sin e ea Meee ot She g Ges oe ee 2 ue Be Be ee Be d ee ee ee ee ee ee ee ee 4 1 5 The case of the single scattering preconditioned integral equation 4 1 6 The case of the Brakhage Werner integral equation 4 1 7 Post processing 2 6 a 4 1 8 Results s s sa 000 ee 4 1 9 Point source wave 2 a a 4 2 The Neumann boundary value problem 02 0004 4 3 Mixing Dirichlet and Neumann boundary conditions 4 4 Penetrable case 2 ee 4 4 1 Integral equation a 0 00000 2 ee ee 4 4 2 A more complex ge
93. rix A Solve the linear system Au b Integral representation Post processing Far field Near field Other incident field displaying geometry Figure 1 Schematic structure of a diff script How to cite p diff Please cite the following reference if you use p diff B Thierry X Antoine C Chniti H Alzubaidi u diff an open Matlab toolbox for computing multiple scattering problems by disks Computer Physics Communications to appear 2015 10 How to install p diff Requirement The toolbox p diff requires the installation of the Matlab software http www mathworks com version 2011 or higher Furthermore j diff should work fine with previous versions of Matlab however without any guarantee Installation The install process is realized as follows 1 Download the p diff toolbox either by using git with git clone http mu diff math cnrs fr git mu diff or by downloading the following zip file and unzipping it in your Matlab working directory or in any other directory that you fix yourself http mu diff math cnrs fr mu diff Download_files mudiff zip Note that git should be prefered to stay up to date easily 2 Add the p diff toolbox to the Matlab s path including subfolders by using either the graphical interface of Matlab or the addpath and save path functions 3 Test the u diff install by typing in the Matlab command window BenchmarkDirichlet This shou
94. s The associated preconditioned EFIE becomes Lp L u p 1 31 22 where the operator L L has the following matrix form T DZIEN o LIDAM 2 2 1 72 1 2 2 172 M z L 2 L I we LAIL LL i i 1 32 LMM 1 M1 LM M 1 M2 p I with I the identity operator Note that this preconditioning accelerates the convergence rate of an iterative solver like e g the GMRES The single scattering preconditioner can be applied to the other integral equations and the following result holds 24 Proposition 1 9 When preconditioned by their single scattering preconditioners the EFIE MFIE and CFIE become identical and similar to the preconditioned BWIE equal up to an invertible operator In other words after being preconditioned the four integral equations lead to the same convergence rate when an iterative solver is applied This proposition shows in particular that there is no need in computing the single scattering preconditioned versions of all the integral equations Because the resulting operator exhibits a good convergence rate for multiple scattering it is hard coded in p diff for both the Dirichlet and the Neumann cases They are based on the single layer potential for the Dirichlet case and the double layer potential for the Neumann case EFIE MFIE or CFIE in both cases 1 3 Mixing Dirichlet and Neumann boundary conditions Let us assume that the collection of M obstacles is such that there is Mp sound soft obstacl
95. s acoustiques et lectromagn tiques Stabilisation d algorithmes it ratifs et aspects de l analyse num rique PhD thesis Universit Paris XI 2006 H Brakhage and P Werner Uber das Dirichletsche Aussenraumproblem fiir die Helmholtzsche Schwingungsgleichung Arch Math 16 325 329 1965 A J Burton and G F Miller The application of integral equation methods to the nu merical solution of some exterior boundary value problems Proc Roy Soc London Ser A 323 201 210 1971 A discussion on numerical analysis of partial differential equations 1970 95 12 13 14 15 16 17 22 23 24 S N Chandler Wilde I G Graham S Langdon and M Lindner Condition number estimates for combined potential boundary integral operators in acoustic scattering J Integral Equations Appl 21 2 229 279 2009 S N Chandler Wilde I G Graham S Langdon and E A Spence Numerical asymptotic boundary integral methods in high frequency acoustic scattering Acta Numerica 21 1 89 305 2012 D L Colton and R Kress Integral Equation Methods in Scattering Theory Pure and Applied Mathematics New York John Wiley amp Sons Inc New York 1983 A Wiley Interscience Publication M Darbas Pr conditionneurs Analytiques de type Calderon pour les Formulations Int grales des Probl mes de Diffraction d Ondes PhD thesis INSA de Toulouse 2004 R Harrington and J Mautz H field E fie
96. tacles atrix_Not_Obstacles MaskMatrixObstacles X Y O a 0 UincOnMesh UincOnMesh Matrix_Not_Obstacles 5 Total field U_tot U UincOnMesh UBWIE_tot UBWIE UincOnMesh Note that to display the near field we suggest to use the following script where the white artefacts are removed thanks to set gcf Renderer Zbuffer and the circles are also displayed using a high zdata here displaying a value called U_tot ind_fig 1 figure ind_fig hold on surf X Y abs U_tot shading interp view 2 colorbar PlotCircles O a ind_fig Color k LineWidth 2 zdata max max abs U_tot set gcf Renderer Zbuffer hold off 4 1 8 Results The figure 4 1 presents the results obtained with the considered configuration that is 3 unit disks placed on 5 2 0 0 and 5 2 with k 1 and an incident plane wave of direction 0 On the figure is shown the obstacles the history of convergence of the GMRES for the 5 integral equations for their dense and sparse storage the radar cross section and also the near fields 4 1 9 Point source wave Our toolbox p diff also provides a right hand side derived from a wave emitted by a point source To solve that kind of problem the right hand side is built thanks to and DnPointSource instead of respectively and DnP laneWave there is currently no equivalent of PlaneWavePrecond For example for the EFIE an
97. tacles must not overlap Useful for example for the points source location a_avoid 1 x N radii of the N holes dmin_avoid 1 x N minimal distance between an obstacle and a hole The holes represented by the _avoid arguments are circular regions where the obstacles must not overlap for example where a point source is emitting a wave Example 1 creating N_scat random disks with random radii is realized by the command O a CreateRandomDisks xmin xmax ymin ymax N_scat amin amax Example 2 building 7 obstacles in the box 10 10 x 10 10 with radii between 0 1 and 0 5 The disks must be separated by a minimal distance equal to 0 1 and without maximal value The command is then O a CreateRandomDisks 10 10 10 10 7 0 1 0 5 0 1 1 Example 3 now consider that a point source is located at 2 2 and that the obstacles must be separated from the source by at least a distance equal to 0 3 Then the _avoid arguments can be used and the resulting function call is O a CreateRandomDisks 10 10 10 10 7 0 1 0 5 0 1 1 2 2 0 3 The disk centered at 2 2 with radius 0 3 is then avoided A second option is to set a_void to zero and set the minimal distance dmin_avoid to 0 3 O a CreateRandomDisks 10 10 10 10 7 0 1 0 5 0 1 1 2 2 0 0 3 Remark 3 1 To check if a disk is correctly placed calls the CheckP l
98. terfaces build the whole vector on only one pattern trace of plane wave normal derivative of a point source wave They are located in the interface directory and their names allow an easy interpretation see also table 3 1 column p diff name DnP laneWave Let us recall that the help informations of the interface functions contains the mathematical description of the incident wave 47 Value y diff function name Param Description 1 PlaneWave beta_inc trace of a plane wave of angle of direction beta_inc A plane wave is defined by e cos 1 sin 8 x2 2 DnPlaneWave beta_inc normal derivative of a plane wave of angle of incidence beta_inc 3 PointSource XS trace of the wave emitted by a point source placed at XS Such a wave is defined in p diff by H klix xsll 4 4 DnPointSource XS normal derivative trace of the wave emitted by a point source placed at Xs 5 PlaneWavePrecond beta_inc same as PlaneWave but multiplied by the inverse of the single layer block diagonal operator see section 2 1 3 on the single scattering preconditioner 6 DnPlaneWavePrecond beta_inc same as DnPlaneWave but multiplied by the inverse of the double layer block diagonal operator Table 3 1 Right hand sides already coded in p diff The two main functions BlockIncidentWave and IncidentWave are now detailed BlockIncidentWave This function computes the block vector of the opposite of the coef
99. the diff user guide license and credits can be found under the directory mudiff Doc 3 2 Common argument and notations A large number of arguments of the y diff functions are similar For the sake of conciseness and if nothing is specified the following arguments refer to the ones given in the Tables below where the indices p and q vary from 1 to M or N_scat in the p diff scripts In addition any function or value specific to u diff is written with the following font Geometry Name Size Description N_scat 1x1 number of obstacles M o 2 x M matrix of the centers of the disks such that 0 1 p resp 0 2 p is the abscissa resp ordinate of the p obstacle Op 2x1 coordinates of the p scatterer Og 2x1 coordinates of the g scatterer a 1 x M vector of the radii of the disks such that a p is the radius of the pt scatterer ap 1x 1 radius of the p scatterer aq 1x1 radius of the g scatterer 40 Common Oc Examples li IntOperators Sparse Block Ful PostProcessing Incident Wave NearField O 2 fe O Preprocessing Incident Wave Figure 3 1 Arborescence of u diff toolbox Parameters wavenumbers incident waves Fourier series expansions Name Size Description beta_inc 1 x1 angle of direction 6 of a plane wave etk cos x1
100. the radar cross section for the different boundary integral equations and the history of convergence of the GMRES solver If you want to see some other examples by using p diff you can try and launch the following time reversal experiments e DORT_NotPenetrable e DORT_dielectric 12 Chapter 1 Boundary integral equations a short survey Contents 1 1 Standard integral equation formulations in acoustic scattering 14 ee ee ee ee ee ee ee ee 14 ed OOM ee eee Ee mH 14 PS gh ORB ee oy OA Radiesse tn Mote ae ae a 17 1 1 4 Brakhage Werner indirect integral equation 19 iy aig Gt aseds gee SEE SRG Ge Ge oS 20 11 6 Summary 42 6 24 4 Se Bei a ee eae Fae ee ee ee 21 sion te christ Oto dor Bb Gace b s Die sets Gets Gone GO Wee es 21 1 2 1 A more explicit writing of the integral equation formulations 21 1 2 2 Single scattering preconditioning 22 aetna ved 23 Oe ere ee eee ee ee ere ee eee ee 24 1 4 1 The boundary value problem 0 2 000 24 sSpeeaes 25 ieee ee be Pes See ete ee wt oe 4 26 The p diff toolbox is based on integral equations and uses the four standard boundary integral operators The derivation and the choice of the integral formulation is let to the user even if some of them are given below In this chapter we provide all the necessary mathematical background to solve time harmonic wave scattering problems by penetrable or impenetrable circular cylinders based on integral equations T
101. ther possibility is to choose A q the interior normal derivative trace Using the traces formulas 1 5 the integral equation 1 14 becomes 1 3 N p nu r 1 16 This Fredholm second kind integral equation named Magnetic Field Integral Equation MFIE is well posed and equivalent to the scattering problem 1 2 as far as k is not an irregular Neumann frequency Proposition 1 6 If k Fn Q7 then the quantity Zp u is the solution of the scattering problem if and only if p is the solution of the MFIE 7 16 For any irregular frequency k of F y Q7 the operator 31 N is no longer one to one In that case and unlike the EFIE the single layer potential p u based on a solution J of the MFIE is not guaranteed to be the solution of the scattering problem 1 2 CFIE Combined Field Integral Equation To avoid the irregular frequencies problem Burton and Miller II considered a linear combi nation of the EFIE and the MFIE by imposing a Fourier Robin boundary condition to up on the boundary T A l a y a 1 17 18 with O0 lt a lt l and Sin 0 1 18 where S 7 is the imaginary part of the complex number 7 Hence the boundary integral equation 1 14 reads as 1 ja a 5 N ont p a a Ogu p onu p 1 19 This Combined Field Integral Equation CFIE denomination of Harrington and Mautz 16 in electromagnetism or Burton Miller integral equation LI is well posed fo
102. twardly directed to Q The domain of propagation is denoted by Q R Q7 When illuminated by an incident time harmonic wave ui the obstacle QT generates a scattered wave field u that is solution to the boundary value problem Aut k u 0 in QT u yine on 1 1 u is outgoing to Q7 The time dependent form of the wave is assumed to be e the wavenumber k 27 w is supposed to be real and positive The first equation of system is the well known Helmholtz equation For the sake of conciseness we consider a Dirichlet boundary condition on I The Neumann case will be studied later The condition u is outgoing to Q means that u satisfies the Sommerfeld s radiation condition at infinity lim x vu ca iku 0 x 00 Ixl where x 27 x2 is the euclidian norm of R Since ui is a solution of the Helmholtz equation in R then the total field up u ui is solution of the following problem Aur keur 0 in QT ur 0 on 1 2 ur ui is outgoing to Q7 It is known that in a suitable mathematical framework the Dirichlet scattering boundary value problems is well posed 14 More precisely we have Theorem 1 1 The problems and 1 2 are uniquely solvable 1 1 2 Volume and boundary integral operators Let the volume single layer integral operator be defined by see e g 20 Theorem 6 12 Z HUD HLR p Zp WER Lox Gyl 14 and the volume double layer integral operato
103. umn with an ordinate equal to 2 5 O a RemoveDisk O_old a_old X 1 Y 2 5 Example 2 Remove the obstacles centered at 2 5 and 3 4 O a RemoveDisk O_old a_old XY 2 3 5 4 Example 3 Create a periodic placement of 11 x 11 unit disks separated by a distance equal to 1 then remove the middle line and the middle column as shown on Figure 2 3 The central disk is moreover centered at 0 0 bx 3 by 3 Nx 11 Ny 11 O RectangularLattice bx by Nx Ny Centered 0 0 a ones l size O 2 O a RemoveDisk O a X 0 Y 0 3 3 2 Truncation of the Fourier series To help the user the formula 2 16 has been coded in p diff The values of Np are stored in a Matlab vector of size 1 x N_scat called M_modes in p diff Computing M_modes only involves the radii of the disks and the wavenumber k which is assumed to be created by the user 46 O0000 OOO0D0O 00000 O0O00O O0000 O0O00O0 OOOO0O0 OO0OO0O0O O0000 O0O00O O0000 O0O00O 0000 O0O000 00000 OOO0N0O OOO00 OO0N0O OOO0O0 OO0O0O 8 ny Figure 3 3 Periodic placement with a row and a column deleted M_modes FourierTruncation a k The resulting vector is such that M_modes p Np where N satisfies 2 16 with obviously a minimal value equal to 0 which consists in only one mode If k is a vector which corresponds to one wavenumber per obstacle
104. value of the scattered field S CFIE 10 BW 0 09 2f Sparse EFIE i Sparse MFIE Sparse CFIE 0 08 Sparse BW 6 val 5 Precond eae Sparse Precond a E 0 06 S 2 k i 8 6f 0 05 o gt 0 Fa 2 0 04 8L P 0 03 6 0 02 104 8 0 01 10 o 12 fi fi fi fi j 10 8 6 4 2 0 2 4 6 8 10 0 5 10 15 20 25 x Iteration number c History of convergence of the GMRES d Absolute value of the scattered field Absolute value of the total field Real part of the total field 10 0 4 0 45 8 0 4 r 0 3 0 35 4 03 K 0 2 0 25 6 on 0 2 2 O15 4 o CA 6 0 05 8 04 o 10 o 8 6 4 2 0 2 4 6 8 10 10 8 6 4 2 0 2 4 6 8 10 x x e Absolute value of the total field f Real part of the total field 10 8 6 4 2 gt 0 2 4 6 8 10 4 Figure 4 2 Different results for a Dirichlet problem for a wave emitted by a point source solved using p diff The first two pictures show respectively the obstacles and the radar cross section while the third one shows the GMRES history of convergence The absolute value of the scattered field is presented on subfigure d and the total field on figures e absolute value and f real part 72 4 2 The Neumann boundary value problem Let us now consider the sound hard scattering problem A k u 0 in QF nu onl u outgoing An efficient solution to this problem is given for example by a preconditioned integral equa
105. where Weight is optional ListOfOperators is a Matlab cell not an array Make use of instead of containing all the sparse matrices that the user wants to involve in the computational process The linear combination of the operators is done at that time for each matrix vector product Example 1 Computing Y LX is done as follows SpL SpIntegralOperator O a M modes k 2 Y SpMatVec X M_modes Spl Example 2 Calculating Y 0 5 x I N X where I is the identity Beware the consists in the following sequence of function calls SpI SpIdentity O a M modes SpN SpDnSingleLayer O a M_modes k Y SpMatVec X M_modes SpI SpN 0 5 1 This can also be done thanks to the SpAddIdentity function as SpN SpDnSingleLayer 0O a M_modes k A SpAddIdentity SpN 0 5 M modes Y SpMatVec X M_modes A Solving a linear system Now that the matrices and the right hand side have been built the sparse linear system can be solved iteratively If one uses for example the GMRES solver 22 the syntax is then Example 1 for the Dirichlet EFIE 1 15 Lp u p SpL SpSingleLayer O a M_modes k Uinc PlaneWave 0O a M modes k beta_inc rho gmres X SpMatVec X M_modes SpL Uinc Example 2 for the Dirichlet MFIE 1 16 1 2 N p nu r SpN SpDnSingleLayer O a M_modes k SpI SpIdentity O a M modes D
106. x I Vu x Vw x dx 1 4 where w q7 and thus satisfies yw vy Since the quantities involved in a scat tering problem do not belong to H Q but rather to A Q the exterior trace and normal derivative trace operators are naturally extended as yf Hb Q gt H L and yE HE Qt A gt H V2 L by yi v af vv and yf v yf vv where v is an arbi trarily compactly supported and infinitely differentiable function on Q which is equal to 1 in a neighborhood of T and where Hj Q A v Hi Q Av L2 Q Remark that when the function v is sufficiently smooth then its normal derivative trace y v given by 1 4 belongs to L P and can be written as yf v x lim eqt_ Vu z n x for almost every x on T Note also that the single and double layer potentials introduced previously belong not only to AL Q U A1 Q but also to HL Q A U HHO A see e g 2 2 Some well known properties of the single and double layer potentials are summarized in the two fol lowing propositions Their proof can be found for example in Theorems 7 5 and 9 6 for proposition and in 20 Theorem 6 12 for proposition 1 3 Proposition 1 2 For any densities p H 2 L and H 2 L the single layer potential Lp and double layer potential MX are outgoing solutions to the Helmholtz equation in R T Moreover the scattered field u solution to can be written as Vx Ent u x
107. xample 2 the Dirichlet BWIE 1 21 the scattered field u is given by u Z Myy and the corresponding far field by eta i k F FarField O a M modes k theta psi eta 1 Example 3 if the densities p and A are known and if the user needs to compute for example M 5 PLpPp PMyXp p 1 then this can be done by the function call F FarField O a M modes k theta rho lambda 1 p 1l p Remark 3 3 Let us remark that the function is also interfaced with the ready to use function FarFieldSingleLayer and FarFieldDoubleLayer These two functions are located in the PostProcessing FarField Interface directory Radar Cross Section RCS The Radar Cross Section RCS can be computed either from the far field with FarField_to_RCs or directly from the density with Rcs In p diff the RCS is obtained by F being the far field o 10log 9 27 F The computation is realized by the function call R FarField_to_RCS F where F has been computed by the FarField function Otherwise the RCS can be computed directly from the densities by R RCS O a M_modes k theta Density Weight 58 where the arguments are exactly the same as for the far field In fact RCS first calls and then Remark 3 4 The RCS function is also interfaced with the two following ready to use func tions RCSSingleLayer and RCSDoubleLayer Both of them ar
108. y IntOperators Sparse Interface Full Sparse matrix of the identity operator SpintegzalOperaton IntOperators Sparse Generic integral operator sparse matrix full SpMat Vec IntOperators Sparse Functions Sparse function sparse matrix vector product possibly multiples SpPrecondDirichlet IntOperators Sparse Interface Full Sparse matrix of the preconditioned integral operator sound soft case SpPrecondNeumann IntOperators Sparse Interface Full Sparse matrix of the preconditioned integral operator sound hard case SpSingleLayer IntOperators Sparse Interface Full Sparse matrix of the single layer integral operator SpSingleMat Vec IntOperators Sparse Functions Sparse function sparse matrix only one vector product e TimeReversalOperator Ezxamples TimeReversal FarField Common Time reversal matrix in acoustic and far field context e TriangularLattice PreProcessing Geometry Build a triangular lattice of disks 87 88 Appendix B List of the u diff functions ordering by folder name Common e repeat_horiz Copy paste a row vector to build a matrix e fangle Angle with horizontal axis e dbesselj First derivative of Bessel function e dbessely First derivative of Newton function e dbesselh First derivative of first kind Hankel function e repeat_vert Copy paste a column vector to build a matrix Examples Benchmark e Example of solution of
109. y V0 0 27 RCS 9 10 logo 27 ax 8 a 7 8 4B To optimize the far fields computation these relations can be written thanks to the inner product between two infinite vectors Indeed let us introduce ay av icp lt m and ay a z i lt p lt m Where ay and av are given by Vp 1 M in 4 a a a 2 E NR e tbpk cos 0 a im 0 r m EVA az P 2 ag co ee i 2 12 Then we obtain the following av ay p and ag 0 av X 2 2 Finite dimensional approximations and numerical solutions proposed in py diff We now have all the ingredients to numerically solve the four integral equations EFIE MFIE CFIE and BWIE for sound soft obstacles In fact any integral equation for any boundary condition can be solved according to the previous developments In practice the resulting infinite Fourier system needs to be truncated to get a finite dimensional problem we must pass 35 from a sum over m E Z to a finite number of Fourier modes that depends on kap p 1 M Let us consider e g the EFIE the extension to the other boundary integral operators being direct The EFIE is given by equation 2 1 Lp U To truncate each Fourier series associated with 5 mez for the obstacle Q we only keep 2N 1 modes in such a way that the indices m of the truncated series satisfy Vp 1 M Np lt m lt Np The truncation parameter N must be fixed large enough with Np gt kap

The µ-diff toolbox: user guide - mu-diff

Contents

Download Pdf Manuals

Related Search

Related Contents

The &micro;-diff toolbox: user guide - mu-diff

Contents

Download Pdf Manuals

Related Search

Related Contents

The µ-diff toolbox: user guide - mu-diff