Implementation in MATLAB of the Iso-Geometric Boundary Elements Method for the resolution of 2D anisotropic elastostatic problems
Abstract
In this work, the implementation of the Isogeometric Boundary Element Method (IGABEM) in MATLAB for the resolution of bi-dimensional anisotropic elastostatic problems is presented. This builds on the degree’s final project of this author, which studied the IGABEM for the isotropic case Firstly, some basics about anisotropic elasticity are given, mostly to clarify nomenclature. Secondly, an introduction to the traditional Boundary Element Method (BEM) is offered, explaining the physics behind it and later particularizing the equations for its numerical implementation. The most important details, or those that take relevance when compared to the IGABEM, are commented. Subsequently, the same process is repeated with the IGABEM, adapting all the equations to this method’s features. The implemented MATLAB code is briefly commented, its difficulties and limitations highlighted. Finally, a comparison between the two methods is done. Analysis are carried out in two elastic problems, one orthotropic and one purely anisotropic, and the results are tested. This does not only show the difference in accuracy but also serves to spot better the differences between these two methods and the advantages of one over the other. Finally, some options for further work are discussed. Together with the ideas, an outline of how they would be made and why they would be useful or interesting is made.
Full text
Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Master’s Degree Final Project Industrial Engineering Implementation in MATLAB of the Iso-Geometric Boundary Elements Method for the resolution of 2D anisotropic elastostatic problems Author: Alberto Sánchez-Reyes González Tutor: Andrés Sáez Pérez Continuum Mechanics and Theory of Structures School of Engineering Universidad de Sevilla Seville, 2020
Master’s Degree Final Project Industrial Engineering Implementation in MATLAB of the Iso-Geometric Boundary Elements Method for the resolution of 2D anisotropic elastostatic problems Author: Alberto Sánchez-Reyes González Tutors: Andrés Sáez Pérez Full Professor Continuum Mechanics and Theory of Structures School of Engineering Universidad de Sevilla Seville, 2020
Master’s Degree Final Project: Implementation in MATLAB of the Iso-Geometric Boundary Elements Method for the resolution of 2D anisotropic elastostatic problems Author: Alberto Sánchez-Reyes González Tutors: Andrés Sáez Pérez El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha:
Abstract In this work, the implementation of the Isogeometric Boundary Element Method (IGABEM) in MATLAB for the resolution of bi-dimensional anisotropic elastostatic problems is presented. This builds on the degree’s final project of this author, which studied the IGABEM for the isotropic case Firstly, some basics about anisotropic elasticity are given, mostly to clarify nomenclature. Secondly, an introduction to the traditional Boundary Element Method (BEM) is offered, explaining the physics behind it and later particularizing the equations for its numerical implementation. The most important details, or those that take relevance when compared to the IGABEM, are commented. Subsequently, the same process is repeated with the IGABEM, adapting all the equations to this method’s features. The implemented MATLAB code is briefly commented, its difficulties and limitations highlighted. Finally, a comparison between the two methods is done. Analysis are carried out in two elastic problems, one orthotropic and one purely anisotropic, and the results are tested. This does not only show the difference in accuracy but also serves to spot better the differences between these two methods and the advantages of one over the other. Finally, some options for further work are discussed. Together with the ideas, an outline of how they would be made and why they would be useful or interesting is made. I
Table of Contents Abstract I List of Figures 1 List of Codes 3 1. Introduction 5 1.1. Motivation 5 1.2. Document organization 6 2. Anisotropic Elasticity 9 2.1. Hooke’s Law 9 2.2. Plane and antiplane problems 10 3. Boundary Integral Equations 11 3.1. Somigliana’s identity 11 3.2. Fundamental solutions 11 3.2.1. Isotropic materials as a degenerate case 12 3.3. Boundary Integral Equations 13 3.3.1. Boundary Points 13 3.3.2. Internal Points 14 Isotropic materials as a degenerate case 14 4. Boundary Element Method 15 4.1. Boundary Elements Discretization 15 4.2. Numerical implementation 17 4.2.1. Rigid body considerations 17 4.2.2. Singular integrals in matrix H 18 4.2.3. Traction discontinuities at corner points 18 4.2.4. Singular integrals in matrix G 18 4.2.5. Code organization 20 ELQUABEANI 20 INPUTEQANI 20 GHMATEQANI 20 EXTINEQANI 21 LOCINEQANI 21 INTEREQANI 21 SIGMAEQANI 22 5. Isogeometric Boundary Element Method 23 5.1. Non-Uniform Rational B-Splines 23 5.1.1. B-Splines 23 III
6Chapter 1. Introduction Figure 1.1 Increase in composite material use in aircraft over past years. Figure taken from [10]. At the cost of having to evaluate more complex functions, several advantages are gained. First and foremost, as has already been stated, the role of the mesh is done by the CAD program in the designing of the part. Besides, this way allows the analysis to be carried out while exactly preserving the geometry of the problem. This takes great importance in problems where the accuracy of the boundary is critically sensitive to geometrical imperfections, such as shell buckling or sliding contact between solids (see Figure 1.3). However, one of the biggest shortcomings of this method is related to the refinement of the analysis. As will be shown in 6.1, sometimes the geometry will be simple (and able to be represented with a small number of functions) while the unknowns fields will not. In such cases, a refinement must be implemented. Nonetheless, the way in which NURBS are built does not allow local refinement, enforcing a global, expensive refinement. An improvement in this field would be the implementation of T-Splines (a superset of NURBS) or PHT-Splines which allow local refinement, solving this deficiency and expanding the method’s adaptivity. For further information, [2] covers T-Splines more deeply and [21] does the same for PHT-splines. This work is centered on two MATLAB codes the author has developed for solving elastic bidimensional problems applying both methods (BEM and IGABEM) to anisotropic materials. The greater part of the document focused on explaining the theory that allow said methods to work; the rest covers the structure of the code and the numerical results it has yielded, which are compared. 1.2 Document organization Chapter 1 is this very section, which makes a case for the value of this work and offers an explained index of its contents. In Chapter 2 some of the ruling equations in anisotropic elasticity will be presented. It will not be an extensive explanation and it will be limited to the equations that are needed to be known in order to understand
1.2 Document organization 7 Figure 1.2 Estimation of the relative time costs of each component of the model generation and analysis process at Sandia National Laboratories. Figure taken from [3]. Figure 1.3 Buckling load of a cylindrical shell with geometrical imperfections. Figure taken from [3]. the following chapters. In Chapter 3 the theoretical explanation behind the equations used in the BEM will be fully developed. The only requisite to the fully understanding of the chapter is to be familiar with the elasticity constitutive equations stated in the prior chapter. In Chapter 4, the traditional BEM will be discussed. First an overall view of how the system of equations is obtained and solved will be presented, then some more specific details about the integrations will be given. For the sake of completeness, the calculation of the fundamental solutions and the free terms will
8Chapter 1. Introduction be included (despite the latter being optional for the execution of the method.) Additionally there will be some comments on how isotropic materials should be treated, being a particular case of the bigger group of anisotropic materials. It will finally introduce the MATLAB program the author has developed without deepening into details. A flux diagram will be attached to show the functioning of the code. In Chapter 5 the IGABEM will be examined. This chapter will start with an introduction to how B-SPLINES and NURBS work. Then the modification of the elasticity equations to fit the IGABEM will be explained. The main problems that arises when doing this will also be commented. In the end, the MATLAB code, analogous to the previous one, will be resumed with its main functionalities presented. Again, a flux diagram is key to the overall comprehension. In Chapter 6 a benchmark example will be shown to compare the accuracy of the different methods. It is a pipe under internal pressure to serve as validation. Next to it a different problem, this time a thin plate with a hole, is solved to test the program in different conditions. Finally, a summary of the presented work is offered, summing up the whole process as well as a brief overall view of the carried out procedure. It is accompanied by a list of possible future works that would improve and extend the validity of this method.
2 Anisotropic Elasticity Since this work will study the behaviour of anisotropic materials, whose properties are usually less known than the isotropic ones, it is convenient to begin by stating the most important equations that rule this behaviour. This topic is discussed with greater detail by F.París, J. Cañas, J.C Marín and A.Barroso [ 6 ], or by T.C.T. Ting in [20]. 2.1 Hooke’s Law The most generic relationship that can be established for a lineal elastic material in 3-D states that: σi j =Ci jkl εkl i,j,k,l=1,2,3(2.1) where repetition of an index implies a summation in that index, σ is the stress tensor, ε means strain and C is the constitutive or elasticity tensor. Even though equation (2.1) establishes 81 different elastic constants, the restrictions about symmetry in the stress tensor can reduce them to 21, for the following properties must be checked: Ci jkl =Cjikl =Ci jlk =Ckli j (2.2) In order to simplify the notation, the following subscripts will be used: i j →αwhereα=iif i=j α=9−i−jif i6=j(2.3) With this in mind the constitutive relations (2.1) can be represented as a matrix: σ11 σ22 σ33 σ23 σ13 σ12 = C11 C12 C13 C14 C15 C16 C22 C23 C24 C25 C26 C33 C34 C35 C36 sym. C44 C45 C46 C55 C56 C66 ε11 ε22 ε33 γ23 γ13 γ12 (2.4) where γi j =2εi j. The inverse relationship to equation 2.4 will be referred as: ε11 ε22 ε33 γ23 γ13 γ12 = S11 S12 S13 S14 S15 S16 S22 S23 S24 S25 S26 S33 S34 S35 S36 sym. S44 S45 S46 S55 S56 S66 σ11 σ22 σ33 σ23 σ13 σ12 (2.5) where Sis known as the compliance tensor. 9
10 Chapter 2. Anisotropic Elasticity If the material has one plane of elastic symmetry, it is called monoclinic and equation 2.4 can be written with 13 constants: σ11 σ22 σ33 σ23 σ13 σ12 = C11 C12 C13 0 0 C16 C22 C23 0 0 C26 C33 0 0 C36 sym. C44 C45 0 C55 0 C66 ε11 ε22 ε33 γ23 γ13 γ12 (2.6) In 2.6 the plane x1−x2has been set as the plane of symmetry. If the material has 3 planes of elastic symmetry, it is called orthotropic and its number of independent constants can be further reduced to 9: σ11 σ22 σ33 σ23 σ13 σ12 = C11 C12 C13 000 C22 C23 000 C33 000 sym. C44 0 0 C55 0 C66 ε11 ε22 ε33 γ23 γ13 γ12 (2.7) Finally, in an isotropic material (where all planes are of symmetry) only 2 constants will remain: σ11 σ22 σ33 σ23 σ13 σ12 = C11 C12 C12 000 C11 C12 000 C11 000 sym. C11−C12 20 0 C11−C12 20 C11−C12 2 ε11 ε22 ε33 γ23 γ13 γ12 (2.8) 2.2 Plane and antiplane problems In plane elasticity, if the material is monoclinic (or one of the cases in equations 2.7 and 2.8), the problem can be split in two separate sets of equations known as the plane and antiplane problems, that can be solved separately. The plane problem is: σ11 σ22 σ12 = C11 C12 C16 C22 C26 sym. C66 ε11 ε22 γ12 (2.9) With ε33 =0 σ33 =C13ε11 +C23ε22 +C36γ12 =for plane strain problems (2.10) ε33 =−C13ε11+C23ε22+C36γ12 C33 σ33 =0)=for plane stress problems (2.11) As for the antiplane problem: σ23 σ13=C44 C45 C45 C55γ23 γ13(2.12)
3 Boundary Integral Equations This chapter will cover the theory behind both BEM and IGABEM methods, presenting the equations that will be particularized and used later. 3.1 Somigliana’s identity The starting point for the mathematical proof of this identity as detailed on the work of Brebbia and Dominguez in [ 5 ]) will be the equilibrium equations of elasticity for a body of domain Ω and boundary σ . Since this work, for simplicity, disregards body forces, such expression is: σk j,j=0in Ω(3.1) In order to minimize (3.1), it can be weighted by a displacement type function u∗ k. Therefore: ZΩ σk j,ju∗ kdΩ=0(3.2) By carrying out an integration by parts twice in the first term of equation (3.2), it is modified to: ZΩ σ∗ k j,jukdΩ+ZΓ pku∗ kdΓ=ZΓ p∗ kukdΓ(3.3) where σ∗ k j and p∗ k are the stresses and tensions related to u∗ k . Now, if the introduced u∗ k function corresponds to the displacement created by and unit force in direction l applied at point y (this function is called fundamental solution), the domain integral only is non null in that particular point, that is: ZΩ σ∗ k j,jukdΩ=ZΩ σ∗ l j,juldΩ=uy l(3.4) As stated before, u∗ and p∗ are called fundamental solutions and will be explained in more detail in the next chapter. The important point here is that among all the arbitrary functions that could have been chosen, u∗ is picked so that the integral in 3.4 becomes a punctual value and p∗ steams from u∗ . Now (3.4) can be substituted into (3.3) by taking each of the separate components of the displacement, and the Somigliana’s identity is reached: uy l+ZΓ p∗ lkukdΓ=ZΓ u∗ lk pkdΓ(3.5) This equation is the basis for the boundary element method. It allows us to work out the displacement in any point of the domain provided the tractions and displacements at the boundaries are known. It is true for any point of the solid, be it internal or on the boundary itself. 3.2 Fundamental solutions In order to apply the above equation (3.5), it is necessary to know the expressions of p∗ lk and u∗ lk . If the x1x2 -plane is taken as the plane of study, p∗ lk and u∗ lk are the tractions and displacements, respectively, created 11
12 Chapter 3. Boundary Integral Equations at point x ( x1,x2 ) when a constant unit load along the x3 -axis in the k direction is applied at point y ( y1,y2 ). The x point is called observation point and y is referred as the collocation point. For an anisotropic bidimensional solid in plane strain and following the nomenclature used by García-Sánchez in [ 7 ], these expressions are (where lis the direction of the component of the solution): u∗ lk =−1 πRe(AkmQmlln(zm−z0 m)(3.6) p∗ lk =1 πRe(LkmQml µmn1−n2 zm−z0 m )(3.7) Where: • The indices l , k , m =1,2 if only the plane problem is being solved (when the material properties allow it, see section 2.2). Otherwise, l,k,m=1,2,3. •Re(x) stands for the real part of x. •nkis the kcomponent of the exterior normal to the boundary Γat point x. •zmis the transformation to the complex plane of the observation point: zm=x1+µmx2(3.8) •z0 mis the transformation to the complex plane of the collocation point: z0 m=y1+µmy2(3.9) •µmare the complex roots with positive imaginary part of the equation: det[C2i j2µ2 m+(C1i j2+C2i j1)µm+C1i j1] = 0(3.10) • The matrix A is made, by columns, from the eigenvectors resulting from solving the characteristic equation of the material: [C2i j2µ2 m+(C1i j2+C2i j1)µm+C1i j1]Ajm =0(no sum in m) (3.11) where the Cαi jβsub-matrices follow the expressions in 2.3. •The matrix Qis defined as: Q=A−1(B−1+B−1)−1(3.12) •The Bmatrix is defined as: B=iAL−1(3.13) •¯ Bstands for the complex conjugate matrix of B. •The matrix Lis made, by columns: Lim = (µmC2i j2+C2i j1)Ajm =−1 µm (C1i j1+µmC1i j2)Ajm (no sum in m) (3.14) It is worth noting that expressions 3.10 to 3.14 depend entirely on the properties of the material, not the geometry, and are constants and not variables. 3.2.1 Isotropic materials as a degenerate case T.C.T Ting explores ([ 20 ]) how the constants in equations 3.10 to 3.14 are depending on the material. It is interesting to point out that µm cannot be real, and that they will come in conjugate pairs. In the particular case of an isotropic material, all roots will be equal and will take the value of i . This will cause the columns of the B matrix to be linearly dependent and its inversion in 3.12 problematic. In the MATLAB program used in this work, the determinant of B is, albeit very small, not exactly zero, which allows it to work for
3.3 Boundary Integral Equations 13 isotropic materials despite being a poorly conditioned problem. The results are acceptably accurate despite this, but it would be advisable to use the fundamental solution for the isotropic case for greater accuracy. For the sake of completeness, they will be covered in this section. The expression of the fundamental solutions (particularized for a 2D isotropic plane strain problem) are: u∗ lk =1 8πµ(1−ν)[(3−4ν)ln 1 rδlk +r,lr,k](3.15) p∗ lk =−1 4π(1−ν)r[r,n((1−2ν)δlk +2r,lr,k)+(1−2ν)(nlr,k−nkr,l)] l,k=1,2(3.16) where r is the distance between the collocation and observation points, µ is the Shear Modulus, ν is the Poisson’s ratio, δlk is the Kronecker delta and r,nis the derivative of rwith respect to the normal, that is: r,n=∂r ∂xk nk+∂r ∂xl nl(3.17) 3.3 Boundary Integral Equations 3.3.1 Boundary Points As expressed in the first section, it is first needed to solve the problem in the boundary to solve the whole domain. Somigliana’s identity is valid for any collocation point, including the boundary. A problem arises, nonetheless, when this happens, since rcauses a singularity when the collocation and integration points match. In order to avoid this singularity, the integration can be carried out around the collocation point, as shown in Figure 3.1. (where Γεis the part of the boundary where the collocation point lies). Figure 3.1 Integration around the singularity. Figure taken from [17]. The equation (3.5) would transform into: uy l+lim ε→0ZΓ−Γε p∗ lkukdΓ+lim ε→0ZΓε p∗ lkuldΓ=lim ε→0ZΓ−Γε u∗ lk pkdΓ+lim ε→0ZΓε u∗ lk pkdΓ(3.18) Only the third and fifth terms bring further trouble. The second and fourth are integrations in the sense of Cauchy principal value. In the fifth, as the limit is approached, pk will become its punctual value and can be placed out of the integral. lim ε→0pkZΓε u∗ lkdΓ(3.19) u∗ lk will produce a term of order ln(1/ε) whereas dΓ is a term of order ε , therefore the value of the limit tends to 0. On the third one the fundamental equation is this time of order ε , the same as the length differential, so
14 Chapter 3. Boundary Integral Equations a non zero coefficient will appear. lim ε→0ukZΓε p∗ lkdΓ=c∗ lkuk(3.20) This term c∗ lk will depend on the shape of the surface where the integration is being carried out as well as the properties of the material. The calculation of these values may prove difficult (especially in 3D) but, as will be seen, it can be avoided. When c∗ lk is moved to the other side of the equation and summed to the displacement term that is there, the result is clk and is typically referred as free term coefficient or jump term coefficient. In conclusion, the equation (3.5) has been transformed into: cy lkuy k+ZΓ p∗ lkukdΓ=ZΓ u∗ lk pkdΓ(3.21) In this expression the integrals are in the sense of Cauchy principal value. This equation will be the central axis of this work and the key to solve firstly the boundary, and secondly the domain. 3.3.2 Internal Points When the collocation point lays inside the solid, the singularity no longer takes place. Therefore equation (3.21) can be used with clk =1 which equals the original (3.5). Known the boundary values, the displacement of any internal point can be easily found. For the computation of stresses, equation 3.21 can be differentiated and the cinematic relationship applied, arriving to an integral equation that allows the computation of stresses: σy sl +ZΓ s∗ slkukdΓ=ZΓ d∗ slk pkdΓ(3.22) where: s∗ slk =Csl jm p∗ jk,m(3.23) d∗ slk =Csl jmu∗ jk,m(3.24) The indices s,l,k,mfollow the same rule as in section 3.2. Isotropic materials as a degenerate case If the same differentiation is made taking the isotropic fundamental solution, the equation equivalent to 3.22 can be written with a slightly different nomenclature, since the derivatives of the fundamental solutions, being simpler, are already included: σy i j +ZΓ s∗ ki jukdΓ=ZΓ d∗ ki j pkdΓ(3.25) where d∗ ki j and s∗ ki j are: d∗ ki j =1 4π(1−ν)r(1−2ν)(δki jr,j+δk jr,i−δi jr,k)+2r,ir,jr,k(3.26) s∗ ki j =µ 2π(1−ν)r22∂r ∂n[(1−2ν)δi jr,k+ν(δikr,j+δjkr,i)−4r,ir,jr,k]+ +2ν(nir,jr,k+njr,ir,k)+(1−2ν)(2nkr,ir,j+njδik +niδjk)−(1−4ν)nkδi j(3.27)
4 Boundary Element Method In this chapter, the traditional BEM will be briefly discussed. It will show how the equations are set and the most important details about their resolution. When the discretization is discussed, it will be generally assumed that only the plane problem is being solved for simplicity regarding indices and sizes. For the most generically anisotropic case, specifics will be mentioned when needed. 4.1 Boundary Elements Discretization In order to solve the boundary, some assumptions must be made. It is going to be split in a number NE of elements with a quantity n of nodes. Since quadratic elements will be used in this work, each element will be formed by three nodes, one in its middle and one in each end as can be seen in Figure 4.1 . Therefore, N=2∗NE. Figure 4.1 Discretization with quadratic elements. Figure taken from [5]. The position, displacement and tractions of any point along an element will be presumed as a combination of the nodal values. This way, shape functions φ(ξ)are introduced, so that: x(ξ) = φ1(ξ)x1+φ2(ξ)x2+φ3(ξ)x3(4.1) u(ξ) = φ1(ξ)u1+φ2(ξ)u2+φ3(ξ)u3(4.2) p(ξ) = φ1(ξ)p1+φ2(ξ)p2+φ3(ξ)p3,(4.3) where xi , ui and pi are the geometry, displacements and tractions, respectively, at node i . ξ is a parametric coordinate which goes over the element varying from -1 to 1, this is, the first node will be placed at ξ=−1 , 15
22 Chapter 4. Boundary Element Method EXTINEQANI for the new integrals (no singularities arise this time since the collocation point is inside the solid, not in the boundary) and SIGMAEQANI. The number of integrations that are carried out are one per internal point per element, which is the number of times the loop in the flow diagram is repeated. After the last integration, the program can move on. SIGMAEQANI This routine gives back the needed terms shown in 3.3.2. These integrals are also computed by means of a standard Gauss quadrature.
5 Isogeometric Boundary Element Method In this chapter the differences and similarities between BEM and IGABEM will be explored. The initial problems raised at working with the IGABEM will be presented and their solution detailed. As introduced in the abstract, NURBS are a special group of functions used by CAD programs to build the boundary of items being designed. The IGABEM has the goal of using these functions in substitution of the shape functions in 4.6 as well as replacing the nodes with a set of points that CAD software programs use. This would ease the step from CAD design to analysis. This apparently simple change has some implications that will be discussed in along the chapter. 5.1 Non-Uniform Rational B-Splines Since NURBS are built from B-Splines, it makes quite sense to introduce them before any step ahead is given. 5.1.1 B-Splines B-Splines are parametric functions that map the parameter space into the physical space. In this work η will be used for the parametric variable. This variable is somehow similar to the ξ used in traditional BEM to define the shape functions. The main difference is that ηis not local to the elements. Possibly the most important components of B-Splines are the knot vector, non-decreasing sequence of coordinates usually referred by Ξ= (η1,η2,...ηn+p+1) where p is the grade the polynomial order (p=2 from now on, for quadratic elements will be employed), n the number of basis functions used and ηi , the knots. This vector splits η into elements. Knot vectors are uniform if the knot are equally distributed in the parameter space, and non-uniform otherwise. They are open if the initial and final knots are repeated p+1 times. This means the splines will be interpolatory at the beginning and end of the curve, and is the most frequent kind of knot vector. The repetition of knot values entails the decrease of continuity order at that point in the spline. In fact, Na_pis p−ltimes differentiable at a knot of multiplicity l. Besides, it is worth noting that the knot vector is indifferent to its scale, that is, the vector (0,0,0,2,2,4,4,4) with η∈[0,4) is equivalent to the vector (0,0,0,1,1,2,2,2) with η∈[0,2) . With this, B-Splines for grade p are defined as: Na_p(η) = η−ηa ηa+p−ηa Na_p−1(η)+ ηa+p+1−η ηa+p+1−ηa+1 Na+1_p−1(η)(5.1) With the exception of grade 0 functions Na_0(η) = 1if ηa≤η<ηa+1 0otherwise (5.2) If equation (5.1) yields the quotient 0 0, it is taken as zero by definition. Figure 5.1 shows an example of B-splines that will be useful for highlighting some facts. 23
24 Chapter 5. Isogeometric Boundary Element Method Figure 5.1 B-splines for the knot vector Ξ= (0,1,2,3,4...). Figure taken from [3]. The first property to be noticed is that, since B-splines are a partition of the unity: n ∑ k=1 Na_p(η) = 1; ∀η(5.3) Secondly, in any point, at most p+1 functions are different from zero. Finally, Na_p(η) = 0 out of the span [ηa,ηa+p+1). 5.1.2 Building NURBS NURBS are built from B-Splines in the following way: Ra_p(η) = Na_p(η)wa ∑n i=1Ni_p(η)wi ,(5.4) where wa is the weight of the function a . It is immediately seen that, if all weights are unit, the NURBS will coincide with the B-splines because of equation (5.3) . Moreover, NURBS also fulfil such equation. The step from B-Splines to NURBS is significant, for it allows to describe certain curves, such as ellipses or circumferences, in an exact way, which was impossible with polynomial functions. With the help of this functions, the coordinates of the curve C described by ηare defined as: C(η) = n ∑ a=1 Ra_p(η)Pa,(5.5) where P are the coordinates of the control points. These control points are similar to the traditional BEM nodes but, very importantly, they do not necessarily belong to the curve. If the first and last control points are coincident, the curve is closed. Figure 5.2 presents an example of a curve built from its control points. The knot vector associated to it is Ξ=0,0,0,1 9,1 9,2 9,1 3,4 9,5 9,2 3,7 3,8 9,8 9,1,1,1 5.2 Boundary Element Discretization In a manner similar to 4.1, the next approximations are obtained:
5.2 Boundary Element Discretization 25 Figure 5.2 Geometry defined by NURBS. Figure taken from [13]. x(η) = n ∑ a=1 Ra_p(η)Pa(5.6) u(η) = n ∑ a=1 Ra_p(η)da(5.7) p(η) = n ∑ a=1 Ra_p(η)qa(5.8) The reason behind having expanded this relations is to draw attention to the d and q vectors. This are called the displacement and traction coefficients, respectively. The term coefficient is used because they do not represent the nodal points displacement or tractions even if such points were all in the boundary (it is worth repeating that the control points do not necessary belong to the curve so they can be placed inside the solid, outside of it or in the boundary). Furthermore, these vectors have no physical meaning. The role of elements will be taken by the different knot spans, provided their length is longer that zero. Analogously to in traditional BEM, the Somigliana’s identity (3.21) can be made discreet as: clk(y) n ∑ a=1 Ra_pdy j+ ne ∑ e=1 n ∑ a=1hZ1 −1 p∗ lk(y,ξ)Ra_p(ξ)|J|(ξ)|J2|dξidae k= = ne ∑ e=1 n ∑ a=1hZ1 −1 u∗ lk(y,ξ)Ra_p(ξ)|J|(ξ)|J2|dξiqae ki,j=1,2(5.9) Here ξ is the same variable it was in traditional BEM; local to the element and ranging from -1 to 1. |J| is also unchanged and defined as in (4.8). |J2| is the Jacobian of the change from η to ξ , since R had been defined respecting η and is used in relation to ξ in this equation. As in 4.7, dae k is the displacement coefficient at control point a of the element e in direction k . In the summations, not all the n elements must be computed, since only p+1will be different from zero. This situation is very similar to the one reached in traditional BEM but now the variables are d and q . After setting the relationship ntimes in different points, the system of equations will be built: H0d=G0q,(5.10) and, after some extra considerations, solved.
26 Chapter 5. Isogeometric Boundary Element Method 5.2.1 Collocation points The first difficulty to face is the choice of the collocation points, that is, where to apply the forces for the fundamental solutions to work. While in traditional BEM the element nodes were the perfect candidates, here their akin control points are not so fitting since they will not generally be in the boundary. Although the control points could be used even though the drawbacks, in this work the Greville abcissae points have been chosen [9], whose coordinates in parametric space are: η0 a=ηa+1+ηa+2+...+ηa+p pa=1,2,...N.(5.11) These points are projections over the boundary of the control ones, as depicted in Figure 5.3: (a) Control points. (b) Collocation points. Figure 5.3 Greville abcissae collocation points. Figure taken from [19]. Again, this is not the only possible option, but it has been chosen following Simpson’s work ([19]). 5.2.2 Boundary conditions The current systems of equations has d and q as unknowns variables. Unlike the previous chapter, where the matrices could be reordered attending to what was given by the boundary conditions, in this case both vectors are completely unknown, for the gathered data is related to u and p . However, (5.7) and (5.8) allow the variables to be transformed. Two matrices A and B can be built, whose terms will be the different NURBS evaluated at the various collocation points, so that: u=Cd (5.12) p=Dq (5.13) which, replaced into (5.10): H0C−1u=G0D−1p(5.14) By joining terms and analogous to (4.9): Hu =Gp (5.15) which can solved in the same way. Both C and D have inverse, since the relation between collocation and control points exists both ways and is unique. 5.2.3 Internal points The computation of internal points was detailed in Equation 3.22. In the IGABEM the integration along the boundary will grant dand qterms, therefore the following adaptation must be made: σy sl +ZΓ s∗ slkCsubdkdΓ=ZΓ d∗ slkDsubqkdΓ(5.16)
5.3 Numerical implementation 27 where Csub are submatrices of C , related to the element or knot span on which the integration is being carried. In this case there would be no difference between Csub and Dsub their differences lies in sharp corners, a single element cannot contain any of them. 5.3 Numerical implementation The numerical implementation of the IGABEM is, in general, very similar to the BEM one. Although the fact that the control and collocation points can lay anywhere makes it harder to handle them, the basic ideas remain the same. In this section, deeper attention will be focused on their differences rather than their similarities. 5.3.1 Rigid body considerations As seen in the previous chapter, the matrix H allowed the computation of its diagonal to be carried out without using the clk terms. However, it is not that clear whether that reasoning can be done in equation (5.14), for d and qhave no physical meaning. Despite this, the rigid body displacement can be applied into (5.14). Since matrix C is made, by rows, of the NURB functions evaluated at the different collocation points, it must sum 1 by rows at stated by 5.3. It can be proven (see Apendix A) than this causes C−1 to sum 1 by rows too. Additionally, it can also be proven (see Apendix A) that this causes H0C−1to sum one by rows. Looking at 5.14 and 5.15, this makes it mandatory for H to sum 1 by rows, and possible for the rigid body consideration to be applied. 5.3.2 Free term coefficient Although the computation of the free term coefficient clk is not mandatory for the method, it has been included for comparing results. The computation of these terms is detailed in [ 16 ], [ 14 ], [ 15 ] and will be included in this work for completeness. As they state in their work, it is limited to materials with a plane of elastic symmetry coincident with the coordinate plane x3=0, in-plane and out-of plane elastic solutions being therefore uncoupled. Only in-plane solutions corresponding to plane strain or plane stress states are analysed here. The equations to be used depend on the nature of the material. If the solid is made of a mathematically non-degenerate material (µmare distinct): ckl =Re 2 ∑ α=1 1 iπκα FlαEkαlogz(1) α z(2) α!(5.17) where: z(e) α=r(e) 1+µαr(e) 2(5.18) E=ξ1(µ1)ξ1(µ2) µ−1 1ξ2(µ1)µ−1 2ξ2(µ2)(5.19) F=−µ1−µ2 1 1 (5.20) κα=2 2 ∑ k=1 EkαFkα(5.21) ξa(µk) = µ2 kSa1−µkSa6+Sa2(5.22) and r(e) i is the i component of the tangent to the boundary on the collocation point, either before the point (e=1) or after it (e=2). The terms Sijare the components of the Smatrix as defined in 2.5. On the contrary, if in case of a mathematically degenerate material (µmare repeated and referred simply as µ) ckl =Re 2 ∑ α,β=1 FlαEkβGαβ (r(e))! 1 2 (5.23)
28 Chapter 5. Isogeometric Boundary Element Method where: E=ξ1(µ1)2η1(µ,¯ µ)−ξ1(µ) µξ1(µ1)¯ µ(2η1(µ,¯ µ)−ξ1(µ))(5.24) ηa(p,q) = pqSa1−0.5(p+q)Sa6+Sa2(5.25) F=−µ−¯ µ 1 1 (5.26) G(x) =1 iπκ ¯z(x)z−1(x)log(z(x)) log(z(x)) 0(5.27) z(x) = x1+µx2 ¯z(x) = x1+¯ µx2(5.28) The other terms (r,S,ξ) are as defined in the non-degenerate case. 5.3.3 Singularities in matrix H The reasoning behind 4.2.2 is no longer valid, since NURBS are not interpolatory at the nodes. In this work, the subtraction of singularity method (SST) has been employed ([ 8 ]). It is a method to carry out singular integrals in the principal Cauchy value sense, and its idea is based on weakening the singularity and splitting the integral into a singular and a regular one. Further details are given in the article. The original article covered the isotropic case for the BEM method, but it has been adapted for the IGABEM anisotropic case. Alternatively, García-Sánchez proposes a different method for dealing with these singularities in [ 7 ] by a change of variable. 5.3.4 Refinement As already mentioned, the NURBS may need an additional refinement to bring valid results. A perfect example is 6.1; in that problem NURBS can represent the geometry of the solid exactly with the use of only 4 elements. However, the displacement and tractions follow a more complex evolution, so a refinement must be made if better accuracy is desired. In this work only h-refinement is treated. P-refinement is a valid option but, since the aim of the work is to compare quadratic examples of both methods, has been not contemplated. This refinement is achieved by inserting new knots in the knot vector. This can be done without altering the curve at all. If the original knot vector was Ξ= (η1,η2,...ηn+p+1) , the new one, called extended knot vector, will be Ξ= (η1=η1,η2,...ηn+m+p+1=ηn+p+1) . The new m+n control points will be linear combinations from the original ones. The way to calculate this new control points and its due weights so the geometry remains unchanged can be consulted in [3]. An example of knot-insertion is depicted in Figure 5.4. 5.3.5 Code organization Now a brief description of the programs will be made. A flux diagram will be adjunct as well, to help the understanding of their behaviour. ELQUABEIGABEMANI This is the main program, which calls the others in order as seen in Figure 5.5. As in the BEM code, it uses cell arrays to separate the data from one boundary to another. Figure 5.5 has two significant differences with the flow diagram that was shown in the BEM chapter (Figure 4.6). Firstly, whenever a singular integral needs to be computed, two different subroutines (LOCINEQANI and SST) will be called. This is because the rigid body motion technique in IGABEM allows to skip the computation of the free terms but not of the singularities in the H0matrix. Secondly, there is an HREFINEMENT subroutine that refines the input boundary curve if needed. This makes the refinement of the mesh much easier than in the BEM case. The rest of the diagram flow is pretty much analogous.
5.3 Numerical implementation 29 Figure 5.4 Knot insertion. Figure taken from [3]. INPUTEQIGABEMANI This subprogram reads a .txt file with the input data. Additionally to what its BEM counterpart did, it also requires a refinement valor (which takes the zero value if no refinement is desired), the knot vector and the weights. Here the disparity between the number of displacement and tractions variable forces the requirement of separate vectors, KODE and KODE1. KODE is related to tractions and takes vale 0 if the traction of that position is unknown, and 1 otherwise. KODE1 is completely analogous, giving information about the displacement.
30 Chapter 5. Isogeometric Boundary Element Method Figure 5.5 Flow diagram IGABEM. HREFINEMENTANI If the refinement value read in INPUTEQIGABEMANI was greater than one, this program inserts knots between the existing ones (refinement+1 knots will take the place of each previous knots). It also computes the new control points and weights to leave the geometry untouched, as commented in 5.3.4. This code is an
5.3 Numerical implementation 31 adaptation of the one provided by [18] BOUNDARY If the control points have been modified, the boundary conditions read by INPUTEQIGABEMANI is no longer enough, for now these new points require their data to be imposed. This subprogram creates the new boundary conditions and is only valid to the particular problem it was written for. GHMATEQIGABEMANI This program computes the different integrations needed for the matrices H0 and G0 . For each of them, calls EXTINEQIGABEMANI if no singularities appear or both SST and LOCINEQIGABEM if they do. After this, computes the matrices C and D , operates with them to reach equation (5.15) and interchanges H and G columns to create the final system. The number of integrations that are carried out are one per collocation point per element, which is the number of times the loop in the flow diagram is repeated. After the last integration, the program can move on. JUMPTERM This program computes the clk terms to be included in the H0matrix. EXTINEQIGABEMANI This program computes the non-singular terms of both matrices using standard Gauss quadrature of ten points. SST This one computes the H0 singular terms using the SST method. The method has been adjusted to tolerate freedom of singularity collocation, since the original document only treats BEM integration. LOCINEQIGABEMANI This program computes the G0 singular terms in the very same way LOCINEQANI did, with a special logarithmic Gauss quadrature. The main difference here that now there are more than three possible node collocation. This has influence in the way the integral is split, but the procedure is similar and not worth repeating. INTEREQIGABEMANI Very similar to its BEM counterpart, this subprogram computes the results of the internal points with the help of EXTINEQIGABEMANI and SIGMAEQIGABEMANI. The number of integrations that are carried out are one per internal point per element, which is the number of times the loop in the flow diagram is repeated. After the last integration, the program can move on. SIGMAEQIGABEMANI Also completely analogous to BEM SIGMAEQANI, this program computes the necessary terms D and S from equations in 3.3.2
38 Chapter 6. Numerical Examples (a) Stress at A. (b) Displacement at A. 6.2 Plate with hole The second problem to be solved is a thin plate with a hole in it. The geometry and loads, symmetry conditions included (minus the movement restriction in the direction perpendicular to the paper, which has not been included in the drawn), are represented in Figure 6.7 where: r=25 mm h =200 mm c =100 mm w=50 mm P =1MPa (6.2)
6.2 Plate with hole 39 (c) Stress at C. (d) Displacement at C. Figure 6.6 Evolution of the results as the mesh is refined.. The material is the same as in the previous example, whose stiffness matrix is (filling out the components that were not specified before), if fully shown (in MPa): C= 138540 3060 0 0 0 0 14550 0 0 0 0 14550 0 0 0 sym. 5840 0 0 5840 0 5840 (6.3) In order to use this material as the most generic anisotropic, it has been rotated 45◦counter-clockwise in
40 Chapter 6. Numerical Examples P P A B r w h c Figure 6.7 Thin plate with hole. regards to each of the three x,y,z axis. The new behaviour matrix, whose transformation is detailed in [ 1 ], is
6.2 Plate with hole 41 (in MPa): C0= 21166 8592 16620 −11037 −10404 7580 22248 17702 −10787 −12184 8121 44519 −21790 −22555 15678 sym. 21312 15296 −10878 22394 −11261 13367 (6.4) The results in MATLABS will again will be compared with those of ANSYS APDL ("Quad 4 node 182" as element type and 8805 total number of elements), both programs working under the assumption of plane stress. The meshes can be seen in Figure 6.8 and Figure 6.9. Figure 6.8 ANSYS mesh. The main principal stress for the line between A and B is shown in Figure 6.10 and Table 6.5. X axis ANSYS BEM Diff IGABEM Diff (m) (MPa) (MPa) (%) (MPa) (%) 0.02 0.45841 0.5396 17.71 0.5162 12.61 0.12 1,7402 1.7587 1.06 1.7059 1.97 0.23 4.1261 4.165 0.94 4.2446 2.87 Table 6.5 1st Principal Stress along the AB line. The total displacement in the plane for the line between A and B is shown in Figure 6.11 and Table 6.6. X axis ANSYS BEM Diff IGABEM Diff (m) (10E-4 m) (10E-4 m) (%) (10E-4 m) (%) 0.02 0.7307 0.7346 0.5364 0.7227 1.0876 0.12 0.7391 0.7408 0.2299 0.7287 1.4020 0.23 0.7678 0.7422 3.3302 0.7356 4.1847 Table 6.6 Total in-plane displacement along the AB line.
42 Chapter 6. Numerical Examples -0.5 0 0.5 1 X axis (dm) 0 0.5 1 1.5 2 Y axis (dm) Boundary Collocation Points Internal Points Figure 6.9 MATLAB mesh. It is interesting to note that, since the program for computing internal points is not prepared to deal with singular integrals, the results become inaccurate when they are too close to the boundary. The distance from which they become unreliable is related to the distance between nodes in the close elements, so a finer mesh allows for smaller precision near the boundary. Figure 6.11 and Table 6.6 show that the accuracy of the solution in displacements decays as the points approach the boundary, although the error doesn’t go over 5%. A refinement in MATLAB has been proven to improve these differences (and the shape of the curves become more like the ANSYS one) although the results being shown have been selected to depict than even a not very dense mesh gives very accurate results away from the boundary and reasonably accurate results close to the boundary. In any way, this example proves that the results of the program are quite good for problems with the most generic anisotropic materials, and with it both plane stress and plane strain have been covered.
6.2 Plate with hole 43 Figure 6.10 1st Principal Stress along the AB line. Figure 6.11 Total in-plane displacement along the AB line.
7 Conclusion and future work In this work, the IGABEM method has been successfully implemented into a MATLAB code to solve practical problems of elasticity. The theory behind the BEM equations have been reasoned, using the elastic equations. The working of the BEM has also been explained with detail, having focused additional attention into the problematic parts of it, such as the singular integrations. Later, these steps have been repeated with the IGABEM. The differences with the BEM have initially presented a difficulty, such as the assignment of collocation points or the application of boundary conditions. These problems have been studied and solved, allowing the IGABEM code to work. There are several interesting results in the IGABEM. The rigid body considerations has been mathematically proved and tested in the code, in combination to the adaptation of the SST method. The features of both methods have been tested. There has not been found a great difference in the accuracy of the solutions, for both methods brought relatively exact results without the need of a large number of elements. The main disparity has been found to lay in the computation and preparation time. While the IGABEM takes longer to compute, its refinement is significantly simpler and easier to do in this program. The second greatest advantage of the IGABEM over the BEM, the precision of the geometry, has not proved to have an decisive impact in the solution of these particular solutions. This is not very surprising since both geometries are quite simple and can be accurately represented with the shape forms in the BEM. In regards to future work, there are several subjects that hold great potential. Firstly, the code could be adapted to work with piezoelectric materials, a subset of anisotropic materials that create an electric voltage when subjected to stress and, inversely, suffer stress when subjected to voltage. Secondly, it could be made to solve 3D problems, so its applicability would be extended. Additionally, there is the possibility to improve the computation of internal points so that they can cope with singular integrals when the points are near the boundary, to avoid loses of accuracy as seen in Section 6.2. Another addition that could improve on the accuracy of this points is the implementation of other splines, such as discussed in the motivation section. This could allow for a local refinement near the internal point, and thus avoid the singularity. 45
Appendix A Demonstrations A.0.1 Property Let C=BA be square matrices of order n. n ∑ j=1 ci j = n ∑ j=1 bi j n ∑ k=1 ajk!(A.1) Demonstration: ci1=bi1a11 +...+binan1 . . . cin =bi1a1n+...+binann (A.2) Summing all terms in A.2 ci1+...+cin =bi1(a11 +...+a1n)+...+bin(an1+...+ann) which proves the property. A.0.2 Conclusion 1 If Aand Bsum 1 by rows, Calso does. Demonstration: n ∑ j=1 ci j = n ∑ j=1 bi j n ∑ k=1 ajk!= n ∑ j=1 bi j =1(A.3) A.0.3 Conclusion 2 If Asums 1 by rows and has an inverse, A−1also does: Demonstration: I=A−1A 1= n ∑ j=1 bi j( n ∑ k=1 ajk) = n ∑ j=1 bi j =1(A.4) 47
54 Chapter B. Codes y{n}(j)=y{n}(j)+Y{n}(k)*R{n}(k,pos); end end figure(NN+1) hold on plot(x{n},y{n}) collocPtsX{n}=zeros(N(n)+1,1); collocPtsY{n}=zeros(N(n)+1,1); collocNURBS=zeros(N(n)+1,N(n)+1); min=1; for k=1:N(n)+1 collocPtsX{n}(k)=0; collocPtsY{n}(k)=0; for j=1:N(n)+1 if j>=min && j<=min+6 collocNURBS(k,j)=NURBSbasis_matlab(j,p,collocPts{n}(k),knotVec{n}, weight{n}); if collocNURBS(k,j)==0 && collocNURBS(k,max(j-1,1))~=0 min=max(j-4,1); break end end collocPtsX{n}(k)= collocPtsX{n}(k)+collocNURBS(k,j)*X{n}(j); collocPtsY{n}(k)= collocPtsY{n}(k)+collocNURBS(k,j)*Y{n}(j); end end hold on plot(collocPtsX{n},collocPtsY{n},'*r') end % Compute the A and B matrices (and also AC); for n=1:NN A{n}=zeros(3*N(n),3*N(n)); AC{n}=zeros(3*N(n),3*N(n)); B{n}=zeros(9*length(active{n}),9*length(active{n})); for i=1:3:3*N(n) for j=1:3:3*N(n) %A{n}(i,j)=NURBSbasis_matlab(ceil(j/2),2,collocPts{n}(ceil(i/2)), knotVec{n},weight{n}); end end for i=2:3:3*N(n) for j=2:3:3*N(n) %A{n}(i,j)=NURBSbasis_matlab(ceil(j/2),2,collocPts{n}(ceil(i/2)), knotVec{n},weight{n}); end end for i=1:3:3*N(n) for j=1:3:3*N(n) %AC{n}(i,j)=NURBSbasis_matlab(ceil(j/2),2,collocPts{n}(ceil(i/2)), knotVec{n},weight{n}); AC{n}(i,j)=collocNURBS(ceil(i/3),ceil(j/3)); AC{n}(i+1,j)=AC{n}(i,j); AC{n}(i+2,j)=AC{n}(i,j); AC{n}(i,j+1)=AC{n}(i,j); AC{n}(i+1,j+1)=AC{n}(i,j);
B.3 GHMATEQIGABEMANI 55 AC{n}(i+2,j+1)=AC{n}(i,j); AC{n}(i,j+2)=AC{n}(i,j); AC{n}(i+1,j+2)=AC{n}(i,j); AC{n}(i+2,j+2)=AC{n}(i,j); A{n}(i,j)=AC{n}(i,j); A{n}(i+1,j+1)=AC{n}(i,j); A{n}(i+2,j+2)=AC{n}(i,j); end end A{n}(end-2,1)=NURBSbasis_matlab(N(n)+1,2,collocPts{n}(end-1),knotVec{n}, weight{n}); A{n}(end-1,2)=NURBSbasis_matlab(N(n)+1,2,collocPts{n}(end-1),knotVec{n}, weight{n}); A{n}(end,3)=NURBSbasis_matlab(N(n)+1,2,collocPts{n}(end-1),knotVec{n}, weight{n}); AC{n}(end-2,1:3)=NURBSbasis_matlab(N(n)+1,2,collocPts{n}(end-1),knotVec{n}, weight{n}); AC{n}(end-1,1:3)=NURBSbasis_matlab(N(n)+1,2,collocPts{n}(end-1),knotVec{n}, weight{n}); AC{n}(end,1:3)=NURBSbasis_matlab(N(n)+1,2,collocPts{n}(end-1),knotVec{n}, weight{n}); for i=1:3:9*length(active{n}) for j=1:3:9*length(active{n}) if ceil(j/3)==activeB{n}(ceil(i/9)) || ceil(j/3)==activeB{n}(ceil(i /9))+1 || ceil(j/3)==activeB{n}(ceil(i/9))+2 B{n}(i,j)=NURBSbasis_matlab(active{n}(ceil(j/9))+floor(rem(j -1,9)/3),2,collocPts{n}(active{n}(ceil(i/9))+floor(rem(i-1,9) /3)),knotVec{n},weight{n}); end end end for i=2:3:9*length(active{n}) for j=2:3:9*length(active{n}) if ceil(j/3)==activeB{n}(ceil(i/9)) || ceil(j/3)==activeB{n}(ceil(i /9))+1 || ceil(j/3)==activeB{n}(ceil(i/9))+2 B{n}(i,j)=NURBSbasis_matlab(active{n}(ceil(j/9))+floor(rem(j -1,9)/3),2,collocPts{n}(active{n}(ceil(i/9))+floor(rem(i-1,9) /3)),knotVec{n},weight{n}); end end end for i=3:3:9*length(active{n}) for j=3:3:9*length(active{n}) if ceil(j/3)==activeB{n}(ceil(i/9)) || ceil(j/3)==activeB{n}(ceil(i /9))+1 || ceil(j/3)==activeB{n}(ceil(i/9))+2 B{n}(i,j)=NURBSbasis_matlab(active{n}(ceil(j/9))+floor(rem(j -1,9)/3),2,collocPts{n}(active{n}(ceil(i/9))+floor(rem(i-1,9) /3)),knotVec{n},weight{n}); end end end if n==1 Afinal=A{n}; ACfinal=AC{n}; Bfinal=B{n}; Activefinal=active{n}; else
56 Chapter B. Codes Afinal=blkdiag(Afinal,A{n}); ACfinal=blkdiag(ACfinal,AC{n}); Bfinal=blkdiag(Bfinal,B{n}); Activefinal=[Activefinal,active{n}+Activefinal(end)+1]; end end G=zeros(3*NT,9*totalactive); H=zeros(3*NT,3*NT); %---------------------------------------------------------------------------------------------------------------% % Compute the GW and HW matrices for each collocation point and each % boundary element. global NURBS dNURBS NURBS=zeros(N(1)+1,(10*size(elRanges{NN},1))); dNURBS=NURBS; pGH=0; for n=1:NN THETA{n}=zeros(N(n),1); for ll=1:N(n) pG=0; pH=0; % Calculate jump term if ll~=1 if ismember(collocPts{n}(ll),uniqueKnots{n}) pos=find(active{n}==ll,1); [C,THETA{n}(ll)]=JUMPTERMANI(X{n}(active{n}(pos)-2),Y{n}(active{n }(pos)-2),X{n}(active{n}(pos)-1),Y{n}(active{n}(pos)-1),X{n}( active{n}(pos)),Y{n}(active{n}(pos)),X{n}(active{n}(pos)),Y{n }(active{n}(pos)),X{n}(active{n}(pos)+1),Y{n}(active{n}(pos) +1),X{n}(active{n}(pos)+2),Y{n}(active{n}(pos)+2),collocPts{n }(ll),collocPts{n}(ll),knotVec{n},weight{n},active{n}(pos)-2, active{n}(pos),mu_k,Ajump,Bjump,K,D); else % C(1,1)=0.5; C(1,2)=0; C(2,1)=0; C(2,2)=0.5; THETA{n}(ll)=pi; end else [C,THETA{n}(ll)]=JUMPTERMANI(X{n}(N(n)-1),Y{n}(N(n)-1),X{n}(N(n)),Y{ n}(N(n)),X{n}(N(n)+1),Y{n}(N(n)+1),X{n}(1),Y{n}(1),X{n}(2),Y{n }(2),X{n}(3),Y{n}(3),collocPts{n}(end),collocPts{n}(1),knotVec{n },weight{n},active{n}(end),active{n}(1),mu_k,Ajump,Bjump,K,D); end % H(2*ll-1+pGH:2*ll+pGH,:)=H(2*ll-1+pGH:2*ll+pGH,:)+repmat(C,[1,NT]) .*ACfinal(2*ll-1+pGH:2*ll+pGH,:); % Carry out the integration for nn=1:NN for i=1:size(elRanges{nn},1) if nn==n if (collocPts{n}(ll)>=elRanges{n}(i,1)&& collocPts{n}(ll)<= elRanges{n}(i,2)) || (ll==1 && i==size(elRanges{nn},1))
B.3 GHMATEQIGABEMANI 57 [HW]=SSTANI(collocPtsX{n}(ll),collocPtsY{n}(ll),X{nn}(active {nn}(i)),Y{nn}(active{nn}(i)),X{nn}(active{nn}(i)+1),Y{ nn}(active{nn}(i)+1),X{nn}(active{nn}(i)+2),Y{nn}(active {nn}(i)+2),elRanges{nn}(i,1),elRanges{nn}(i,2),collocPts {n}(ll),knotVec{nn},weight{nn},active{nn}(i),Qmat,Lmat, mu_k,collocPtsX{n}(1),collocPtsY{n}(1),i); [GW]=LOCINEQIGABEMANI2(collocPtsX{n}(ll),collocPtsY{n}(ll),X {nn}(active{nn}(i)),Y{nn}(active{nn}(i)),X{nn}(active{nn }(i)+1),Y{nn}(active{nn}(i)+1),X{nn}(active{nn}(i)+2),Y{ nn}(active{nn}(i)+2),elRanges{n}(i,1),elRanges{n}(i,2), collocPts{n}(ll),knotVec{n},weight{n},active{nn}(i),Qmat ,Amat,mu_k,collocPtsX{n}(1),collocPtsY{n}(1),i); else [HW,GW]=EXTINEQIGABEMANI(collocPtsX{n}(ll),collocPtsY{n}(ll) ,X{nn}(active{nn}(i)),Y{nn}(active{nn}(i)),X{nn}(active{ nn}(i)+1),Y{nn}(active{nn}(i)+1),X{nn}(active{nn}(i)+2), Y{nn}(active{nn}(i)+2),elRanges{nn}(i,1),elRanges{nn}(i ,2),knotVec{nn},weight{nn},active{nn}(i),Qmat,Amat,Lmat, mu_k,collocPtsX{n}(1),collocPtsY{n}(1),i); end else [HW,GW]=EXTINEQIGABEMANI(collocPtsX{n}(ll),collocPtsY{n}(ll),X{ nn}(active{nn}(i)),Y{nn}(active{nn}(i)),X{nn}(active{nn}(i) +1),Y{nn}(active{nn}(i)+1),X{nn}(active{nn}(i)+2),Y{nn}( active{nn}(i)+2),elRanges{nn}(i,1),elRanges{nn}(i,2), knotVec{nn},weight{nn},active{nn}(i),Qmat,Amat,Lmat,mu_k, collocPtsX{n}(1),collocPtsY{n}(1),i); end % Plug the GW and the HW matrices into the general G and H matrices. for k=1:3 for j=1:9 G(3*ll-3+k+pGH,4.5*(active{nn}(i)-1)+j+pG)=G(3*ll-3+k+pGH ,4.5*(active{nn}(i)-1)+j+pG)+GW(k,j); if (i==size(elRanges{nn},1) && j>=7) H(3*ll-3+k+pGH,j+pH-6)=H(3*ll-3+k+pGH,j+pH-6)+HW(k,j); else H(3*ll-3+k+pGH,3*(active{nn}(i)-1)+j+pH)=H(3*ll-3+k+ pGH,3*(active{nn}(i)-1)+j+pH)+HW(k,j); end end end if i==n1{nn}(1) pG=pG+9*size(elRanges{nn},1); pH=pH+3*N(nn); end end end if ll==N(n) pGH=pGH+3*ll; end end end
58 Chapter B. Codes %------------------------------------------------------------------------------------------------------------------------------------------------------------------ % If the jump term is not in use, apply rigid movement: for i=1:1:3*NT sum1=sum(H(i,1:3:end)); sum2=sum(H(i,2:3:end)); sum3=sum(H(i,3:3:end)); for j=1:3:3*NT H(i,j)=H(i,j)-sum1*AC{1}(i,j); end for j=2:3:3*NT H(i,j)=H(i,j)-sum2*AC{1}(i,j); end for j=3:3:3*NT H(i,j)=H(i,j)-sum3*AC{1}(i,j); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Hfinal=H/Afinal; Gfinal=G/Bfinal; sum(Hfinal') % Create KODE and DFI as usable arrays. for k=1:sum(N) KODEu(3*k-2)=displacements(6*k-5); KODEu(3*k-1)=displacements(6*k-3); KODEu(3*k)=displacements(6*k-1); DFIu(3*k-2)=displacements(6*k-4); DFIu(3*k-1)=displacements(6*k-2); DFIu(3*k)=displacements(6*k); end KODEu=logical(KODEu); DFIu=DFIu(KODEu); for k=1:length(tractions)/2 KODEv(k)=tractions(2*k-1); DFIv(k)=tractions(2*k); end % Manage the variables of G (and B) for n=1:NN if n==1 Thetafinal=THETA{n}; else Thetafinal=[Thetafinal,THETA{n}]; end end for j=1:size(Gfinal,2)
B.3 GHMATEQIGABEMANI 59 Gcolumn(j)=Activefinal(ceil(j/9))+floor(rem(j-1,9)/3); end pos=0; for n=1:NN Gcolumn(9*length(active{n})-2+pos:9*length(active{n})+pos)=Gcolumn(9*length (active{n})-2+pos:9*length(active{n})+pos)-3*size(elRanges{nn},1); pos=9*length(active{n})+pos; end Gs=Gcolumn; for j=10:3:size(Gfinal,2) for k=1:3:j-3 if Gcolumn(j)~=0 if Gcolumn(k)==Gcolumn(j) && (pi-20*eps)<=Thetafinal(Gcolumn(j)) && Thetafinal(Gcolumn(j))<=(pi+20*eps) Gfinal(:,k)=Gfinal(:,j)+Gfinal(:,k); Gfinal(:,k+1)=Gfinal(:,j+1)+Gfinal(:,k+1); Gfinal(:,k+2)=Gfinal(:,j+2)+Gfinal(:,k+2); Gfinal(:,j+2)=0; Gfinal(:,j+1)=0; Gfinal(:,j)=0; Gcolumn(:,j+2)=0; Gcolumn(:,j+1)=0; Gcolumn(:,j)=0; end end end end KODE1v=KODEv; DFI1v=DFIv; for j=size(Gcolumn,2):-1:1 if Gfinal(:,j)==0 Gfinal(:,j)=[]; Gcolumn(:,j)=[]; KODEv(j)=[]; DFIv(j)=[]; end end KODEv=logical(KODEv); DFIv=DFIv(KODEv); % Reorder the columns of the system of equations in accordance with the % boundary conditions and form system matrix C*x=FI. Du=Hfinal(:,KODEu); Cu=Hfinal(:,~KODEu); Dv=Gfinal(:,KODEv); Cv=Gfinal(:,~KODEv); D=[-Du,Dv]; C=[Cu,-Cv]; DFI=[DFIu,DFIv];
60 Chapter B. Codes FI=D*DFI'; B.4 JUMPTERMANI function [C,THETA]=JUMPTERMANI(X1,Y1,X2,Y2,X3,Y3,X4,Y4,X5,Y5,X6,Y6,COL1,COL2, knotVec,weight,active1,active2,mu_k,Ajump,Bjump,K,D) [~,dR11]=NURBSbasis_matlab(active1,2,COL1-10*eps,knotVec,weight); [~,dR21]=NURBSbasis_matlab(active1+1,2,COL1-10*eps,knotVec,weight); [~,dR31]=NURBSbasis_matlab(active1+2,2,COL1-10*eps,knotVec,weight); [~,dR12]=NURBSbasis_matlab(active2,2,COL2+10*eps,knotVec,weight); [~,dR22]=NURBSbasis_matlab(active2+1,2,COL2+10*eps,knotVec,weight); [~,dR32]=NURBSbasis_matlab(active2+2,2,COL2+10*eps,knotVec,weight); A1=dR11*X1+dR21*X2+dR31*X3; B1=dR11*Y1+dR21*Y2+dR31*Y3; MG1=sqrt(A1^2+B1^2); R11=-A1/MG1; R21=-B1/MG1; A2=dR12*X4+dR22*X5+dR32*X6; B2=dR12*Y4+dR22*Y5+dR32*Y6; MG2=sqrt(A2^2+B2^2); R12=A2/MG2; R22=B2/MG2; if D==1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% z1=R11+mu_k(1)*R21; z2=R12+mu_k(1)*R22; G1=1/(pi*1i*K)*[conj(z1)/z1,log(z1);log(z1),0]; G2=1/(pi*1i*K)*[conj(z2)/z2,log(z2);log(z2),0]; if COL1==3 G1=1/(pi*1i*K)*[conj(z1)/z1,log(z1)+2*pi*1i;log(z1)+2*pi*1i,0]; end C(1,1)=real(Bjump(1,1)*Ajump(1,1)*G1(1,1)+Bjump(1,1)*Ajump(1,2)*G1(1,2)+ Bjump(1,2)*Ajump(1,1)*G1(2,1)+Bjump(1,2)*Ajump(1,2)*G1(2,2)-Bjump(1,1)* Ajump(1,1)*G2(1,1)-Bjump(1,1)*Ajump(1,2)*G2(1,2)-Bjump(1,2)*Ajump(1,1)* G2(2,1)-Bjump(1,2)*Ajump(1,2)*G2(2,2)); C(1,2)=real(Bjump(1,1)*Ajump(2,1)*G1(1,1)+Bjump(1,1)*Ajump(2,2)*G1(1,2)+ Bjump(1,2)*Ajump(2,1)*G1(2,1)+Bjump(1,2)*Ajump(2,2)*G1(2,2)-Bjump(1,1)* Ajump(2,1)*G2(1,1)-Bjump(1,1)*Ajump(2,2)*G2(1,2)-Bjump(1,2)*Ajump(2,1)* G2(2,1)-Bjump(1,2)*Ajump(2,2)*G2(2,2)); C(2,1)=real(Bjump(2,1)*Ajump(1,1)*G1(1,1)+Bjump(2,1)*Ajump(1,2)*G1(1,2)+ Bjump(2,2)*Ajump(1,1)*G1(2,1)+Bjump(2,2)*Ajump(1,2)*G1(2,2)-Bjump(2,1)* Ajump(1,1)*G2(1,1)-Bjump(2,1)*Ajump(1,2)*G2(1,2)-Bjump(2,2)*Ajump(1,1)* G2(2,1)-Bjump(2,2)*Ajump(1,2)*G2(2,2)); C(2,2)=real(Bjump(2,1)*Ajump(2,1)*G1(1,1)+Bjump(2,1)*Ajump(2,2)*G1(1,2)+ Bjump(2,2)*Ajump(2,1)*G1(2,1)+Bjump(2,2)*Ajump(2,2)*G1(2,2)-Bjump(2,1)*
B.5 LOCINEQIGABEMANI 61 Ajump(2,1)*G2(1,1)-Bjump(2,1)*Ajump(2,2)*G2(1,2)-Bjump(2,2)*Ajump(2,1)* G2(2,1)-Bjump(2,2)*Ajump(2,2)*G2(2,2)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% else z11=R11+mu_k(1)*R21; z12=R12+mu_k(1)*R22; z21=R11+mu_k(2)*R21; z22=R12+mu_k(2)*R22; C(1,1)=real(1/(pi*1i*K(1))*Bjump(1,1)*Ajump(1,1)*log(z11/z12)+1/(pi*1i*K(2) )*Bjump(1,2)*Ajump(1,2)*log(z21/z22)); C(1,2)=real(1/(pi*1i*K(1))*Bjump(1,1)*Ajump(2,1)*log(z11/z12)+1/(pi*1i*K(2) )*Bjump(1,2)*Ajump(2,2)*log(z21/z22)); C(2,1)=real(1/(pi*1i*K(1))*Bjump(2,1)*Ajump(1,1)*log(z11/z12)+1/(pi*1i*K(2) )*Bjump(2,2)*Ajump(1,2)*log(z21/z22)); C(2,2)=real(1/(pi*1i*K(1))*Bjump(2,1)*Ajump(2,1)*log(z11/z12)+1/(pi*1i*K(2) )*Bjump(2,2)*Ajump(2,2)*log(z21/z22)); end THETA=atan2d(A2*-B1-B2*-A1,A2*-A1+B2*-B1)*pi/180; if THETA<0 THETA=THETA+2*pi; end if THETA<0 THETA=THETA+2*pi; end end B.5 LOCINEQIGABEMANI function [GW]=LOCINEQIGABEMANI2(XP,YP,X1,Y1,X2,Y2,X3,Y3,inicio,final,COL, knotVec,weight,active,Qmat,Amat,mu_k,Xfirst,Yfirst,vecpos) %This subroutine computes the GW matriz when the collocation point is one of the nodes of the integration element. %The coefficients are computed by numerical integration: the non singular %part is computed using standard Gauss quadrature, the logarithmic part is %computed using a special quadrature formula. global NURBS dNURBS %Data for the Gauss quadrature. GI=[0.9739065285,-0.9739065285,0.8650633666,-0.8650633666,0.6794095682, -0.6794095682,0.4333953941,-0.4333953941,0.1488743389,-0.1488743389]; OME=[0.0666713443,0.0666713443,0.1494513491,0.1494513491,0.2190863625, 0.2190863625,0.2692667193,0.2692667193,0.2955242247,0.2955242247]; % Data for the special quadrature.
62 Chapter B. Codes GIL=[0.0090426309,0.0539712662,0.1353118246,0.2470524162,0.3802125396, 0.5237923179,0.6657752055,0.7941904160,0.8981610912,0.9688479887]; OMEL=[0.1209551319,0.1863635425,0.1956608732,0.1735771421,0.1356956729, 0.0936467585,0.0557877273,0.0271598109,0.0095151826,0.0016381576]; if COL<inicio COL=final; end if COL~=inicio && COL~=final NODO=2; elseif COL==inicio NODO=1; elseif COL==final NODO=3; end GW=zeros(3,9); JAC=(final-inicio)/2; % Transformation to complex plane zo1=XP+YP*mu_k(1); zo2=XP+YP*mu_k(2); zo3=XP+YP*mu_k(3); for i=1:10 % Compute shape functions for numerical integration. if XP==Xfirst && YP==Yfirst [F1,dR1]=NURBSbasis_matlab(active,2,(GI(i)*(final-inicio)+final+inicio) /2,knotVec,weight); [F2,dR2]=NURBSbasis_matlab(active+1,2,(GI(i)*(final-inicio)+final+ inicio)/2,knotVec,weight); [F3,dR3]=NURBSbasis_matlab(active+2,2,(GI(i)*(final-inicio)+final+ inicio)/2,knotVec,weight); NURBS(active,10*vecpos-10+i)=F1; NURBS(active+1,10*vecpos-10+i)=F2; NURBS(active+2,10*vecpos-10+i)=F3; dNURBS(active,10*vecpos-10+i)=dR1; dNURBS(active+1,10*vecpos-10+i)=dR2; dNURBS(active+2,10*vecpos-10+i)=dR3; else F1=NURBS(active,10*vecpos-10+i); F2=NURBS(active+1,10*vecpos-10+i); F3=NURBS(active+2,10*vecpos-10+i); dR1=dNURBS(active,10*vecpos-10+i); dR2=dNURBS(active+1,10*vecpos-10+i); dR3=dNURBS(active+2,10*vecpos-10+i); end XCO=X1*F1+X2*F2+X3*F3; YCO=Y1*F1+Y2*F2+Y3*F3; % Transformation to complex plane
B.5 LOCINEQIGABEMANI 63 z1=XCO+YCO*mu_k(1); z2=XCO+YCO*mu_k(2); z3=XCO+YCO*mu_k(3); % Jacobian A=dR1*X1+dR2*X2+dR3*X3; B=dR1*Y1+dR2*Y2+dR3*Y3; XJA=sqrt(A^2+B^2); %Compute GW coefficients if NODO==1 [FL11,dR11]=NURBSbasis_matlab(active,2,GIL(i)*(final-inicio)+inicio, knotVec,weight); [FL12,dR12]=NURBSbasis_matlab(active+1,2,GIL(i)*(final-inicio)+inicio, knotVec,weight); [FL13,dR13]=NURBSbasis_matlab(active+2,2,GIL(i)*(final-inicio)+inicio, knotVec,weight); AL=dR11*X1+dR12*X2+dR13*X3; BL=dR11*Y1+dR12*Y2+dR13*Y3; XJA2=sqrt(AL^2+BL^2); L11=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(1,1)*Qmat(1,1)+Amat(1,2)*Qmat(2,1) +Amat(1,3)*Qmat(3,1)); L12=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(2,1)*Qmat(1,1)+Amat(2,2)*Qmat(2,1) +Amat(2,3)*Qmat(3,1)); L13=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(3,1)*Qmat(1,1)+Amat(3,2)*Qmat(2,1) +Amat(3,3)*Qmat(3,1)); L21=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(1,1)*Qmat(1,2)+Amat(1,2)*Qmat(2,2) +Amat(1,3)*Qmat(3,2)); L22=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(2,1)*Qmat(1,2)+Amat(2,2)*Qmat(2,2) +Amat(2,3)*Qmat(3,2)); L23=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(3,1)*Qmat(1,2)+Amat(3,2)*Qmat(2,2) +Amat(3,3)*Qmat(3,2)); L31=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(1,1)*Qmat(1,3)+Amat(1,2)*Qmat(2,3) +Amat(1,3)*Qmat(3,3)); L32=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(2,1)*Qmat(1,3)+Amat(2,2)*Qmat(2,3) +Amat(2,3)*Qmat(3,3)); L33=+JAC*XJA2*2*OMEL(i)/pi*real(Amat(3,1)*Qmat(1,3)+Amat(3,2)*Qmat(2,3) +Amat(3,3)*Qmat(3,3)); NU=(GI(i)+1)/2; S11=-JAC*XJA*OME(i)/pi*real(Amat(1,1)*Qmat(1,1)*log((z1-zo1)/NU)+Amat (1,2)*Qmat(2,1)*log((z2-zo2)/NU)+Amat(1,3)*Qmat(3,1)*log((z3-zo3)/ NU)); S12=-JAC*XJA*OME(i)/pi*real(Amat(2,1)*Qmat(1,1)*log((z1-zo1)/NU)+Amat (2,2)*Qmat(2,1)*log((z2-zo2)/NU)+Amat(2,3)*Qmat(3,1)*log((z3-zo3)/ NU)); S13=-JAC*XJA*OME(i)/pi*real(Amat(3,1)*Qmat(1,1)*log((z1-zo1)/NU)+Amat (3,2)*Qmat(2,1)*log((z2-zo2)/NU)+Amat(3,3)*Qmat(3,1)*log((z3-zo3)/ NU));
70 Chapter B. Codes % HW(1,1)=HW(1,1)+OME(i)*XJA*F1/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(1,2)=HW(1,2)+OME(i)*XJA*F1/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(2,1)=HW(2,1)+OME(i)*XJA*F1/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(2,2)=HW(2,2)+OME(i)*XJA*F1/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % GW(1,3)=GW(1,3)-OME(i)*XJA*F2/pi*real(Amat(1,1)*Qmat(1,1)*log(z1-zo1)+ Amat(1,2)*Qmat(2,1)*log(z2-zo2)); % GW(1,4)=GW(1,4)-OME(i)*XJA*F2/pi*real(Amat(2,1)*Qmat(1,1)*log(z1-zo1)+ Amat(2,2)*Qmat(2,1)*log(z2-zo2)); % GW(2,3)=GW(2,3)-OME(i)*XJA*F2/pi*real(Amat(1,1)*Qmat(1,2)*log(z1-zo1)+ Amat(1,2)*Qmat(2,2)*log(z2-zo2)); % GW(2,4)=GW(2,4)-OME(i)*XJA*F2/pi*real(Amat(2,1)*Qmat(1,2)*log(z1-zo1)+ Amat(2,2)*Qmat(2,2)*log(z2-zo2)); % HW(1,3)=HW(1,3)+OME(i)*XJA*F2/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(1,4)=HW(1,4)+OME(i)*XJA*F2/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(2,3)=HW(2,3)+OME(i)*XJA*F2/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(2,4)=HW(2,4)+OME(i)*XJA*F2/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % GW(1,5)=GW(1,5)-OME(i)*XJA*F3/pi*real(Amat(1,1)*Qmat(1,1)*log(z1-zo1)+ Amat(1,2)*Qmat(2,1)*log(z2-zo2)); % GW(1,6)=GW(1,6)-OME(i)*XJA*F3/pi*real(Amat(2,1)*Qmat(1,1)*log(z1-zo1)+ Amat(2,2)*Qmat(2,1)*log(z2-zo2)); % GW(2,5)=GW(2,5)-OME(i)*XJA*F3/pi*real(Amat(1,1)*Qmat(1,2)*log(z1-zo1)+ Amat(1,2)*Qmat(2,2)*log(z2-zo2)); % GW(2,6)=GW(2,6)-OME(i)*XJA*F3/pi*real(Amat(2,1)*Qmat(1,2)*log(z1-zo1)+ Amat(2,2)*Qmat(2,2)*log(z2-zo2)); % HW(1,5)=HW(1,5)+OME(i)*XJA*F3/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(1,6)=HW(1,6)+OME(i)*XJA*F3/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(2,5)=HW(2,5)+OME(i)*XJA*F3/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); % HW(2,6)=HW(2,6)+OME(i)*XJA*F3/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_k (1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)); GW(1,1)=GW(1,1)-OME(i)*XJA*F1/pi*real(Amat(1,1)*Qmat(1,1)*log(z1-zo1)+Amat (1,2)*Qmat(2,1)*log(z2-zo2)+Amat(1,3)*Qmat(3,1)*log(z3-zo3)); GW(1,2)=GW(1,2)-OME(i)*XJA*F1/pi*real(Amat(2,1)*Qmat(1,1)*log(z1-zo1)+Amat (2,2)*Qmat(2,1)*log(z2-zo2)+Amat(2,3)*Qmat(3,1)*log(z3-zo3)); GW(1,3)=GW(1,3)-OME(i)*XJA*F1/pi*real(Amat(3,1)*Qmat(1,1)*log(z1-zo1)+Amat (3,2)*Qmat(2,1)*log(z2-zo2)+Amat(3,3)*Qmat(3,1)*log(z3-zo3)); GW(1,4)=GW(1,4)-OME(i)*XJA*F2/pi*real(Amat(1,1)*Qmat(1,1)*log(z1-zo1)+Amat (1,2)*Qmat(2,1)*log(z2-zo2)+Amat(1,3)*Qmat(3,1)*log(z3-zo3)); GW(1,5)=GW(1,5)-OME(i)*XJA*F2/pi*real(Amat(2,1)*Qmat(1,1)*log(z1-zo1)+Amat (2,2)*Qmat(2,1)*log(z2-zo2)+Amat(2,3)*Qmat(3,1)*log(z3-zo3)); GW(1,6)=GW(1,6)-OME(i)*XJA*F2/pi*real(Amat(3,1)*Qmat(1,1)*log(z1-zo1)+Amat (3,2)*Qmat(2,1)*log(z2-zo2)+Amat(3,3)*Qmat(3,1)*log(z3-zo3)); GW(1,7)=GW(1,7)-OME(i)*XJA*F3/pi*real(Amat(1,1)*Qmat(1,1)*log(z1-zo1)+Amat (1,2)*Qmat(2,1)*log(z2-zo2)+Amat(1,3)*Qmat(3,1)*log(z3-zo3)); GW(1,8)=GW(1,8)-OME(i)*XJA*F3/pi*real(Amat(2,1)*Qmat(1,1)*log(z1-zo1)+Amat (2,2)*Qmat(2,1)*log(z2-zo2)+Amat(2,3)*Qmat(3,1)*log(z3-zo3));
B.6 EXTINEQIGABEMANI 71 GW(1,9)=GW(1,9)-OME(i)*XJA*F3/pi*real(Amat(3,1)*Qmat(1,1)*log(z1-zo1)+Amat (3,2)*Qmat(2,1)*log(z2-zo2)+Amat(3,3)*Qmat(3,1)*log(z3-zo3)); GW(2,1)=GW(2,1)-OME(i)*XJA*F1/pi*real(Amat(1,1)*Qmat(1,2)*log(z1-zo1)+Amat (1,2)*Qmat(2,2)*log(z2-zo2)+Amat(1,3)*Qmat(3,2)*log(z3-zo3)); GW(2,2)=GW(2,2)-OME(i)*XJA*F1/pi*real(Amat(2,1)*Qmat(1,2)*log(z1-zo1)+Amat (2,2)*Qmat(2,2)*log(z2-zo2)+Amat(2,3)*Qmat(3,2)*log(z3-zo3)); GW(2,3)=GW(2,3)-OME(i)*XJA*F1/pi*real(Amat(3,1)*Qmat(1,2)*log(z1-zo1)+Amat (3,2)*Qmat(2,2)*log(z2-zo2)+Amat(3,3)*Qmat(3,2)*log(z3-zo3)); GW(2,4)=GW(2,4)-OME(i)*XJA*F2/pi*real(Amat(1,1)*Qmat(1,2)*log(z1-zo1)+Amat (1,2)*Qmat(2,2)*log(z2-zo2)+Amat(1,3)*Qmat(3,2)*log(z3-zo3)); GW(2,5)=GW(2,5)-OME(i)*XJA*F2/pi*real(Amat(2,1)*Qmat(1,2)*log(z1-zo1)+Amat (2,2)*Qmat(2,2)*log(z2-zo2)+Amat(2,3)*Qmat(3,2)*log(z3-zo3)); GW(2,6)=GW(2,6)-OME(i)*XJA*F2/pi*real(Amat(3,1)*Qmat(1,2)*log(z1-zo1)+Amat (3,2)*Qmat(2,2)*log(z2-zo2)+Amat(3,3)*Qmat(3,2)*log(z3-zo3)); GW(2,7)=GW(2,7)-OME(i)*XJA*F3/pi*real(Amat(1,1)*Qmat(1,2)*log(z1-zo1)+Amat (1,2)*Qmat(2,2)*log(z2-zo2)+Amat(1,3)*Qmat(3,2)*log(z3-zo3)); GW(2,8)=GW(2,8)-OME(i)*XJA*F3/pi*real(Amat(2,1)*Qmat(1,2)*log(z1-zo1)+Amat (2,2)*Qmat(2,2)*log(z2-zo2)+Amat(2,3)*Qmat(3,2)*log(z3-zo3)); GW(2,9)=GW(2,9)-OME(i)*XJA*F3/pi*real(Amat(3,1)*Qmat(1,2)*log(z1-zo1)+Amat (3,2)*Qmat(2,2)*log(z2-zo2)+Amat(3,3)*Qmat(3,2)*log(z3-zo3)); GW(3,1)=GW(3,1)-OME(i)*XJA*F1/pi*real(Amat(1,1)*Qmat(1,3)*log(z1-zo1)+Amat (1,2)*Qmat(2,3)*log(z2-zo2)+Amat(1,3)*Qmat(3,3)*log(z3-zo3)); GW(3,2)=GW(3,2)-OME(i)*XJA*F1/pi*real(Amat(2,1)*Qmat(1,3)*log(z1-zo1)+Amat (2,2)*Qmat(2,3)*log(z2-zo2)+Amat(2,3)*Qmat(3,3)*log(z3-zo3)); GW(3,3)=GW(3,3)-OME(i)*XJA*F1/pi*real(Amat(3,1)*Qmat(1,3)*log(z1-zo1)+Amat (3,2)*Qmat(2,3)*log(z2-zo2)+Amat(3,3)*Qmat(3,3)*log(z3-zo3)); GW(3,4)=GW(3,4)-OME(i)*XJA*F2/pi*real(Amat(1,1)*Qmat(1,3)*log(z1-zo1)+Amat (1,2)*Qmat(2,3)*log(z2-zo2)+Amat(1,3)*Qmat(3,3)*log(z3-zo3)); GW(3,5)=GW(3,5)-OME(i)*XJA*F2/pi*real(Amat(2,1)*Qmat(1,3)*log(z1-zo1)+Amat (2,2)*Qmat(2,3)*log(z2-zo2)+Amat(2,3)*Qmat(3,3)*log(z3-zo3)); GW(3,6)=GW(3,6)-OME(i)*XJA*F2/pi*real(Amat(3,1)*Qmat(1,3)*log(z1-zo1)+Amat (3,2)*Qmat(2,3)*log(z2-zo2)+Amat(3,3)*Qmat(3,3)*log(z3-zo3)); GW(3,7)=GW(3,7)-OME(i)*XJA*F3/pi*real(Amat(1,1)*Qmat(1,3)*log(z1-zo1)+Amat (1,2)*Qmat(2,3)*log(z2-zo2)+Amat(1,3)*Qmat(3,3)*log(z3-zo3)); GW(3,8)=GW(3,8)-OME(i)*XJA*F3/pi*real(Amat(2,1)*Qmat(1,3)*log(z1-zo1)+Amat (2,2)*Qmat(2,3)*log(z2-zo2)+Amat(2,3)*Qmat(3,3)*log(z3-zo3)); GW(3,9)=GW(3,9)-OME(i)*XJA*F3/pi*real(Amat(3,1)*Qmat(1,3)*log(z1-zo1)+Amat (3,2)*Qmat(2,3)*log(z2-zo2)+Amat(3,3)*Qmat(3,3)*log(z3-zo3)); HW(1,1)=HW(1,1)+OME(i)*XJA*F1/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,2)=HW(1,2)+OME(i)*XJA*F1/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,3)=HW(1,3)+OME(i)*XJA*F1/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,4)=HW(1,4)+OME(i)*XJA*F2/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,5)=HW(1,5)+OME(i)*XJA*F2/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,6)=HW(1,6)+OME(i)*XJA*F2/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2));
72 Chapter B. Codes HW(1,7)=HW(1,7)+OME(i)*XJA*F3/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,8)=HW(1,8)+OME(i)*XJA*F3/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,9)=HW(1,9)+OME(i)*XJA*F3/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,1)=HW(2,1)+OME(i)*XJA*F1/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,2)=HW(2,2)+OME(i)*XJA*F1/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,3)=HW(2,3)+OME(i)*XJA*F1/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,4)=HW(2,4)+OME(i)*XJA*F2/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,5)=HW(2,5)+OME(i)*XJA*F2/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,6)=HW(2,6)+OME(i)*XJA*F2/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,7)=HW(2,7)+OME(i)*XJA*F3/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,8)=HW(2,8)+OME(i)*XJA*F3/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,9)=HW(2,9)+OME(i)*XJA*F3/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,1)=HW(3,1)+OME(i)*XJA*F1/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,2)=HW(3,2)+OME(i)*XJA*F1/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,3)=HW(3,3)+OME(i)*XJA*F1/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,4)=HW(3,4)+OME(i)*XJA*F2/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,5)=HW(3,5)+OME(i)*XJA*F2/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,6)=HW(3,6)+OME(i)*XJA*F2/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,7)=HW(3,7)+OME(i)*XJA*F3/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(1,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2));
B.7 SSTANI 73 HW(3,8)=HW(3,8)+OME(i)*XJA*F3/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(2,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,9)=HW(3,9)+OME(i)*XJA*F3/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)*(mu_k(1) *ETA1-ETA2)+Lmat(3,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)* Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); end GW=GW/JAC; HW=HW/JAC; end B.7 SSTANI function [HW]=SSTANI(XP,YP,X1,Y1,X2,Y2,X3,Y3,inicio,final,COL,knotVec,weight, active,Qmat,Lmat,mu_k,Xfirst,Yfirst,vecpos) if COL<inicio COL=final; end global NURBS dNURBS GI=[0.9739065285,-0.9739065285,0.8650633665,-0.8650633665,0.6794095683, -0.6794095683,0.4333953941,-0.4333953941,0.1488743389,-0.1488743389]; OME=[0.0666713443,0.0666713443,0.1494513491,0.1494513491,0.2190863625, 0.2190863625,0.2692667193,0.2692667193,0.2955242247,0.2955242247]; HW=zeros(3,9); JAC=2/(final-inicio); NODO=2; if COL~=inicio && COL~=final NODO=2; elseif COL==inicio NODO=1; elseif COL==final NODO=3; end % Transformation to complex plane zo1=XP+YP*mu_k(1); zo2=XP+YP*mu_k(2); zo3=XP+YP*mu_k(3); if NODO==1 || NODO==3 for i=1:10 % Compute the values of the shape functions at the integration points. if XP==Xfirst && YP==Yfirst [F1,dR1]=NURBSbasis_matlab(active,2,(GI(i)*(final-inicio)+final+inicio)/2, knotVec,weight);
74 Chapter B. Codes [F2,dR2]=NURBSbasis_matlab(active+1,2,(GI(i)*(final-inicio)+final+inicio) /2,knotVec,weight); [F3,dR3]=NURBSbasis_matlab(active+2,2,(GI(i)*(final-inicio)+final+inicio) /2,knotVec,weight); NURBS(active,10*vecpos-10+i)=F1; NURBS(active+1,10*vecpos-10+i)=F2; NURBS(active+2,10*vecpos-10+i)=F3; dNURBS(active,10*vecpos-10+i)=dR1; dNURBS(active+1,10*vecpos-10+i)=dR2; dNURBS(active+2,10*vecpos-10+i)=dR3; else F1=NURBS(active,10*vecpos-10+i); F2=NURBS(active+1,10*vecpos-10+i); F3=NURBS(active+2,10*vecpos-10+i); dR1=dNURBS(active,10*vecpos-10+i); dR2=dNURBS(active+1,10*vecpos-10+i); dR3=dNURBS(active+2,10*vecpos-10+i); end % Compute geometrical properties at the integration points. XCO=X1*F1+X2*F2+X3*F3; YCO=Y1*F1+Y2*F2+Y3*F3; A=dR1*X1+dR2*X2+dR3*X3; B=dR1*Y1+dR2*Y2+dR3*Y3; XJA=sqrt(A^2+B^2); ETA1=B/XJA; ETA2=-A/XJA; % Transformation to complex plane z1=XCO+YCO*mu_k(1); z2=XCO+YCO*mu_k(2); z3=XCO+YCO*mu_k(3); % Compute GW and HW matrices. HW(1,1)=HW(1,1)+OME(i)*XJA/JAC*F1/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,2)=HW(1,2)+OME(i)*XJA/JAC*F1/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,3)=HW(1,3)+OME(i)*XJA/JAC*F1/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,4)=HW(1,4)+OME(i)*XJA/JAC*F2/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,5)=HW(1,5)+OME(i)*XJA/JAC*F2/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,6)=HW(1,6)+OME(i)*XJA/JAC*F2/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2));
B.7 SSTANI 75 HW(1,7)=HW(1,7)+OME(i)*XJA/JAC*F3/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,8)=HW(1,8)+OME(i)*XJA/JAC*F3/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(1,9)=HW(1,9)+OME(i)*XJA/JAC*F3/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,1)=HW(2,1)+OME(i)*XJA/JAC*F1/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,2)=HW(2,2)+OME(i)*XJA/JAC*F1/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,3)=HW(2,3)+OME(i)*XJA/JAC*F1/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,4)=HW(2,4)+OME(i)*XJA/JAC*F2/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,5)=HW(2,5)+OME(i)*XJA/JAC*F2/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,6)=HW(2,6)+OME(i)*XJA/JAC*F2/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,7)=HW(2,7)+OME(i)*XJA/JAC*F3/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,8)=HW(2,8)+OME(i)*XJA/JAC*F3/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(2,9)=HW(2,9)+OME(i)*XJA/JAC*F3/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,1)=HW(3,1)+OME(i)*XJA/JAC*F1/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,2)=HW(3,2)+OME(i)*XJA/JAC*F1/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,3)=HW(3,3)+OME(i)*XJA/JAC*F1/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,4)=HW(3,4)+OME(i)*XJA/JAC*F2/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,5)=HW(3,5)+OME(i)*XJA/JAC*F2/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,6)=HW(3,6)+OME(i)*XJA/JAC*F2/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,7)=HW(3,7)+OME(i)*XJA/JAC*F3/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2));
76 Chapter B. Codes HW(3,8)=HW(3,8)+OME(i)*XJA/JAC*F3/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); HW(3,9)=HW(3,9)+OME(i)*XJA/JAC*F3/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); end if NODO==1 [F1,dR11]=NURBSbasis_matlab(active,2,inicio+(final-inicio)*0.0001,knotVec, weight); [F2,dR12]=NURBSbasis_matlab(active+1,2,inicio+(final-inicio)*0.0001,knotVec ,weight); [F3,dR13]=NURBSbasis_matlab(active+2,2,inicio+(final-inicio)*0.0001,knotVec ,weight); A=dR11*X1+dR12*X2+dR13*X3; B=dR11*Y1+dR12*Y2+dR13*Y3; XJA=sqrt(A^2+B^2); XCO=X1*F1+X2*F2+X3*F3; YCO=Y1*F1+Y2*F2+Y3*F3; z1=XCO+YCO*mu_k(1); z2=XCO+YCO*mu_k(2); z3=XCO+YCO*mu_k(3); ETA1=B/XJA; ETA2=-A/XJA; CONST11=XJA*(final-inicio)*0.0001/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST12=XJA*(final-inicio)*0.0001/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST13=XJA*(final-inicio)*0.0001/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST21=XJA*(final-inicio)*0.0001/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST22=XJA*(final-inicio)*0.0001/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST23=XJA*(final-inicio)*0.0001/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST31=XJA*(final-inicio)*0.0001/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST32=XJA*(final-inicio)*0.0001/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST33=XJA*(final-inicio)*0.0001/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); end
B.7 SSTANI 77 if NODO==3 [F1,dR31]=NURBSbasis_matlab(active,2,final-(final-inicio)*0.0001,knotVec, weight); [F2,dR32]=NURBSbasis_matlab(active+1,2,final-(final-inicio)*0.0001,knotVec, weight); [F3,dR33]=NURBSbasis_matlab(active+2,2,final-(final-inicio)*0.0001,knotVec, weight); A=dR31*X1+dR32*X2+dR33*X3; B=dR31*Y1+dR32*Y2+dR33*Y3; XJA=sqrt(A^2+B^2); XCO=X1*F1+X2*F2+X3*F3; YCO=Y1*F1+Y2*F2+Y3*F3; z1=XCO+YCO*mu_k(1); z2=XCO+YCO*mu_k(2); ETA1=B/XJA; ETA2=-A/XJA; CONST11=XJA*(final-inicio)*0.0001/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST12=XJA*(final-inicio)*0.0001/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST13=XJA*(final-inicio)*0.0001/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,1)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,1)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST21=XJA*(final-inicio)*0.0001/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST22=XJA*(final-inicio)*0.0001/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST23=XJA*(final-inicio)*0.0001/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,2)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,2)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST31=XJA*(final-inicio)*0.0001/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(1,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (1,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST32=XJA*(final-inicio)*0.0001/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(2,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (2,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2)); CONST33=XJA*(final-inicio)*0.0001/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)*(mu_ k(1)*ETA1-ETA2)+Lmat(3,2)*Qmat(2,3)/(z2-zo2)*(mu_k(2)*ETA1-ETA2)+Lmat (3,3)*Qmat(3,3)/(z3-zo3)*(mu_k(3)*ETA1-ETA2));; end HW(1,1)=HW(1,1)+F1*log(abs(XJA/JAC))*CONST11; HW(1,2)=HW(1,2)+F1*log(abs(XJA/JAC))*CONST12; HW(1,3)=HW(1,3)+F1*log(abs(XJA/JAC))*CONST13; HW(1,4)=HW(1,4)+F2*log(abs(XJA/JAC))*CONST11; HW(1,5)=HW(1,5)+F2*log(abs(XJA/JAC))*CONST12; HW(1,6)=HW(1,6)+F2*log(abs(XJA/JAC))*CONST13; HW(1,7)=HW(1,7)+F3*log(abs(XJA/JAC))*CONST11; HW(1,8)=HW(1,8)+F3*log(abs(XJA/JAC))*CONST12; HW(1,9)=HW(1,9)+F3*log(abs(XJA/JAC))*CONST13;
78 Chapter B. Codes HW(2,1)=HW(2,1)+F1*log(abs(XJA/JAC))*CONST21; HW(2,2)=HW(2,2)+F1*log(abs(XJA/JAC))*CONST22; HW(2,3)=HW(2,3)+F1*log(abs(XJA/JAC))*CONST23; HW(2,4)=HW(2,4)+F2*log(abs(XJA/JAC))*CONST21; HW(2,5)=HW(2,5)+F2*log(abs(XJA/JAC))*CONST22; HW(2,6)=HW(2,6)+F2*log(abs(XJA/JAC))*CONST23; HW(2,7)=HW(2,7)+F3*log(abs(XJA/JAC))*CONST21; HW(2,8)=HW(2,8)+F3*log(abs(XJA/JAC))*CONST22; HW(2,9)=HW(2,9)+F3*log(abs(XJA/JAC))*CONST23; HW(3,1)=HW(3,1)+F1*log(abs(XJA/JAC))*CONST31; HW(3,2)=HW(3,2)+F1*log(abs(XJA/JAC))*CONST32; HW(3,3)=HW(3,3)+F1*log(abs(XJA/JAC))*CONST33; HW(3,4)=HW(3,4)+F2*log(abs(XJA/JAC))*CONST31; HW(3,5)=HW(3,5)+F2*log(abs(XJA/JAC))*CONST32; HW(3,6)=HW(3,6)+F2*log(abs(XJA/JAC))*CONST33; HW(3,7)=HW(3,7)+F3*log(abs(XJA/JAC))*CONST31; HW(3,8)=HW(3,8)+F3*log(abs(XJA/JAC))*CONST32; HW(3,9)=HW(3,9)+F3*log(abs(XJA/JAC))*CONST33; end if NODO==2 for i=1:10 [F11,dR11]=NURBSbasis_matlab(active,2,(GI(i)*(COL-inicio)+COL+inicio)/2, knotVec,weight); [F21,dR21]=NURBSbasis_matlab(active+1,2,(GI(i)*(COL-inicio)+COL+inicio)/2, knotVec,weight); [F31,dR31]=NURBSbasis_matlab(active+2,2,(GI(i)*(COL-inicio)+COL+inicio)/2, knotVec,weight); [F12,dR12]=NURBSbasis_matlab(active,2,(GI(i)*(final-COL)+COL+final)/2, knotVec,weight); [F22,dR22]=NURBSbasis_matlab(active+1,2,(GI(i)*(final-COL)+COL+final)/2, knotVec,weight); [F32,dR32]=NURBSbasis_matlab(active+2,2,(GI(i)*(final-COL)+COL+final)/2, knotVec,weight); XCO1=X1*F11+X2*F21+X3*F31; YCO1=Y1*F11+Y2*F21+Y3*F31; A1=dR11*X1+dR21*X2+dR31*X3; B1=dR11*Y1+dR21*Y2+dR31*Y3; XJA1=sqrt(A1^2+B1^2); ETA11=B1/XJA1; ETA12=-A1/XJA1; XCO2=X1*F12+X2*F22+X3*F32; YCO2=Y1*F12+Y2*F22+Y3*F32; A2=dR12*X1+dR22*X2+dR32*X3; B2=dR12*Y1+dR22*Y2+dR32*Y3; XJA2=sqrt(A2^2+B2^2); ETA21=B2/XJA2; ETA22=-A2/XJA2; z11=XCO1+YCO1*mu_k(1); z12=XCO1+YCO1*mu_k(2); z13=XCO1+YCO1*mu_k(3); z21=XCO2+YCO2*mu_k(1);
B.7 SSTANI 79 z22=XCO2+YCO2*mu_k(2); z23=XCO2+YCO2*mu_k(3); HW(1,1)=HW(1,1)+OME(i)*(COL-inicio)/2*XJA1*F11/pi*real(Lmat(1,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(1,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(1,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F12/pi*real(Lmat(1,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(1,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(1,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22)); HW(1,2)=HW(1,2)+OME(i)*(COL-inicio)/2*XJA1*F11/pi*real(Lmat(2,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(2,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(2,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F12/pi*real(Lmat(2,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(2,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(2,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22)); HW(1,3)=HW(1,3)+OME(i)*(COL-inicio)/2*XJA1*F11/pi*real(Lmat(3,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(3,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(3,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F12/pi*real(Lmat(3,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(3,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(3,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22)); HW(1,4)=HW(1,4)+OME(i)*(COL-inicio)/2*XJA1*F21/pi*real(Lmat(1,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(1,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(1,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F22/pi*real(Lmat(1,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(1,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(1,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22)); HW(1,5)=HW(1,5)+OME(i)*(COL-inicio)/2*XJA1*F21/pi*real(Lmat(2,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(2,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(2,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F22/pi*real(Lmat(2,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(2,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(2,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22)); HW(1,6)=HW(1,6)+OME(i)*(COL-inicio)/2*XJA1*F21/pi*real(Lmat(3,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(3,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(3,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F22/pi*real(Lmat(3,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(3,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(3,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22)); HW(1,7)=HW(1,7)+OME(i)*(COL-inicio)/2*XJA1*F31/pi*real(Lmat(1,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(1,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(1,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F32/pi*real(Lmat(1,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(1,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(1,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22)); HW(1,8)=HW(1,8)+OME(i)*(COL-inicio)/2*XJA1*F31/pi*real(Lmat(2,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(2,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(2,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F32/pi*real(Lmat(2,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(2,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(2,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22)); HW(1,9)=HW(1,9)+OME(i)*(COL-inicio)/2*XJA1*F31/pi*real(Lmat(3,1)*Qmat(1,1) /(z11-zo1)*(mu_k(1)*ETA11-ETA12)+Lmat(3,2)*Qmat(2,1)/(z12-zo2)*(mu_k(2) *ETA11-ETA12)+Lmat(3,3)*Qmat(3,1)/(z13-zo3)*(mu_k(3)*ETA11-ETA12))+OME( i)*(final-COL)/2*XJA2*F32/pi*real(Lmat(3,1)*Qmat(1,1)/(z21-zo1)*(mu_k (1)*ETA21-ETA22)+Lmat(3,2)*Qmat(2,1)/(z22-zo2)*(mu_k(2)*ETA21-ETA22)+ Lmat(3,3)*Qmat(3,1)/(z23-zo3)*(mu_k(3)*ETA21-ETA22));
86 Chapter B. Codes U21_2=1/pi*real(Amat(1,1)*Qmat(1,2)/(z1-zo1)*mu_k(1)+Amat(1,2)*Qmat(2,2)/(z 2-zo2)*mu_k(2)+Amat(1,3)*Qmat(3,2)/(z3-zo3)*mu_k(3)); U22_1=1/pi*real(Amat(2,1)*Qmat(1,2)/(z1-zo1)+Amat(2,2)*Qmat(2,2)/(z2-zo2)+ Amat(2,3)*Qmat(3,2)/(z3-zo3)); U22_2=1/pi*real(Amat(2,1)*Qmat(1,2)/(z1-zo1)*mu_k(1)+Amat(2,2)*Qmat(2,2)/(z 2-zo2)*mu_k(2)+Amat(2,3)*Qmat(3,2)/(z3-zo3)*mu_k(3)); U23_1=1/pi*real(Amat(3,1)*Qmat(1,2)/(z1-zo1)+Amat(3,2)*Qmat(2,2)/(z2-zo2)+ Amat(3,3)*Qmat(3,2)/(z3-zo3)); U23_2=1/pi*real(Amat(3,1)*Qmat(1,2)/(z1-zo1)*mu_k(1)+Amat(3,2)*Qmat(2,2)/(z 2-zo2)*mu_k(2)+Amat(3,3)*Qmat(3,2)/(z3-zo3)*mu_k(3)); U31_1=1/pi*real(Amat(1,1)*Qmat(1,3)/(z1-zo1)+Amat(1,2)*Qmat(2,3)/(z2-zo2)+ Amat(1,3)*Qmat(3,3)/(z3-zo3)); U31_2=1/pi*real(Amat(1,1)*Qmat(1,3)/(z1-zo1)*mu_k(1)+Amat(1,2)*Qmat(2,3)/(z 2-zo2)*mu_k(2)+Amat(1,3)*Qmat(3,3)/(z3-zo3)*mu_k(3)); U32_1=1/pi*real(Amat(2,1)*Qmat(1,3)/(z1-zo1)+Amat(2,2)*Qmat(2,3)/(z2-zo2)+ Amat(2,3)*Qmat(3,3)/(z3-zo3)); U32_2=1/pi*real(Amat(2,1)*Qmat(1,3)/(z1-zo1)*mu_k(1)+Amat(2,2)*Qmat(2,3)/(z 2-zo2)*mu_k(2)+Amat(2,3)*Qmat(3,3)/(z3-zo3)*mu_k(3)); U33_1=1/pi*real(Amat(3,1)*Qmat(1,3)/(z1-zo1)+Amat(3,2)*Qmat(2,3)/(z2-zo2)+ Amat(3,3)*Qmat(3,3)/(z3-zo3)); U33_2=1/pi*real(Amat(3,1)*Qmat(1,3)/(z1-zo1)*mu_k(1)+Amat(3,2)*Qmat(2,3)/(z 2-zo2)*mu_k(2)+Amat(3,3)*Qmat(3,3)/(z3-zo3)*mu_k(3)); P11_1=1/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (1,2)*Qmat(2,1)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)*Qmat(3,1)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2)); P11_2=1/pi*real(Lmat(1,1)*Qmat(1,1)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(1,2)*Qmat(2,1)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(1,3)* Qmat(3,1)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); P12_1=1/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (2,2)*Qmat(2,1)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)*Qmat(3,1)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2)); P12_2=1/pi*real(Lmat(2,1)*Qmat(1,1)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(2,2)*Qmat(2,1)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(2,3)* Qmat(3,1)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); P13_1=1/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (3,2)*Qmat(2,1)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)*Qmat(3,1)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2)); P13_2=1/pi*real(Lmat(3,1)*Qmat(1,1)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(3,2)*Qmat(2,1)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(3,3)* Qmat(3,1)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); P21_1=1/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (1,2)*Qmat(2,2)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)*Qmat(3,2)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2)); P21_2=1/pi*real(Lmat(1,1)*Qmat(1,2)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(1,2)*Qmat(2,2)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(1,3)* Qmat(3,2)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); P22_1=1/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (2,2)*Qmat(2,2)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)*Qmat(3,2)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2)); P22_2=1/pi*real(Lmat(2,1)*Qmat(1,2)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(2,2)*Qmat(2,2)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(2,3)* Qmat(3,2)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); P23_1=1/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (3,2)*Qmat(2,2)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)*Qmat(3,2)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2));
B.9 SIGMAEQIGABEMANI 87 P23_2=1/pi*real(Lmat(3,1)*Qmat(1,2)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(3,2)*Qmat(2,2)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(3,3)* Qmat(3,2)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); P31_1=1/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (1,2)*Qmat(2,3)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(1,3)*Qmat(3,3)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2)); P31_2=1/pi*real(Lmat(1,1)*Qmat(1,3)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(1,2)*Qmat(2,3)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(1,3)* Qmat(3,3)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); P32_1=1/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (2,2)*Qmat(2,3)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(2,3)*Qmat(3,3)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2)); P32_2=1/pi*real(Lmat(2,1)*Qmat(1,3)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(2,2)*Qmat(2,3)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(2,3)* Qmat(3,3)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); P33_1=1/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)+Lmat (3,2)*Qmat(2,3)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)+Lmat(3,3)*Qmat(3,3)/(z3- zo3)^2*(mu_k(3)*ETA1-ETA2)); P33_2=1/pi*real(Lmat(3,1)*Qmat(1,3)/(z1-zo1)^2*(mu_k(1)*ETA1-ETA2)*mu_k(1)+ Lmat(3,2)*Qmat(2,3)/(z2-zo2)^2*(mu_k(2)*ETA1-ETA2)*mu_k(2)+Lmat(3,3)* Qmat(3,3)/(z3-zo3)^2*(mu_k(3)*ETA1-ETA2)*mu_k(3)); S11(1)=S11(1)+OME(i)*XJA*F1*(Cmat(1,1)*P11_1+Cmat(1,6)*P11_2+Cmat(1,6)*P 21_1+Cmat(1,2)*P21_2+Cmat(1,5)*P31_1+Cmat(1,4)*P31_2); S11(2)=S11(2)+OME(i)*XJA*F1*(Cmat(1,1)*P12_1+Cmat(1,6)*P12_2+Cmat(1,6)*P 22_1+Cmat(1,2)*P22_2+Cmat(1,5)*P32_1+Cmat(1,4)*P32_2); S11(3)=S11(3)+OME(i)*XJA*F1*(Cmat(1,1)*P13_1+Cmat(1,6)*P13_2+Cmat(1,6)*P 23_1+Cmat(1,2)*P23_2+Cmat(1,5)*P33_1+Cmat(1,4)*P33_2); S11(4)=S11(4)+OME(i)*XJA*F2*(Cmat(1,1)*P11_1+Cmat(1,6)*P11_2+Cmat(1,6)*P 21_1+Cmat(1,2)*P21_2+Cmat(1,5)*P31_1+Cmat(1,4)*P31_2); S11(5)=S11(5)+OME(i)*XJA*F2*(Cmat(1,1)*P12_1+Cmat(1,6)*P12_2+Cmat(1,6)*P 22_1+Cmat(1,2)*P22_2+Cmat(1,5)*P32_1+Cmat(1,4)*P32_2); S11(6)=S11(6)+OME(i)*XJA*F2*(Cmat(1,1)*P13_1+Cmat(1,6)*P13_2+Cmat(1,6)*P 23_1+Cmat(1,2)*P23_2+Cmat(1,5)*P33_1+Cmat(1,4)*P33_2); S11(7)=S11(7)+OME(i)*XJA*F3*(Cmat(1,1)*P11_1+Cmat(1,6)*P11_2+Cmat(1,6)*P 21_1+Cmat(1,2)*P21_2+Cmat(1,5)*P31_1+Cmat(1,4)*P31_2); S11(8)=S11(8)+OME(i)*XJA*F3*(Cmat(1,1)*P12_1+Cmat(1,6)*P12_2+Cmat(1,6)*P 22_1+Cmat(1,2)*P22_2+Cmat(1,5)*P32_1+Cmat(1,4)*P32_2); S11(9)=S11(9)+OME(i)*XJA*F3*(Cmat(1,1)*P13_1+Cmat(1,6)*P13_2+Cmat(1,6)*P 23_1+Cmat(1,2)*P23_2+Cmat(1,5)*P33_1+Cmat(1,4)*P33_2); S12(1)=S12(1)+OME(i)*XJA*F1*(Cmat(6,1)*P11_1+Cmat(6,6)*P11_2+Cmat(6,6)*P 21_1+Cmat(6,2)*P21_2+Cmat(6,5)*P31_1+Cmat(6,4)*P31_2); S12(2)=S12(2)+OME(i)*XJA*F1*(Cmat(6,1)*P12_1+Cmat(6,6)*P12_2+Cmat(6,6)*P 22_1+Cmat(6,2)*P22_2+Cmat(6,5)*P32_1+Cmat(6,4)*P32_2); S12(3)=S12(3)+OME(i)*XJA*F1*(Cmat(6,1)*P13_1+Cmat(6,6)*P13_2+Cmat(6,6)*P 23_1+Cmat(6,2)*P23_2+Cmat(6,5)*P33_1+Cmat(6,4)*P33_2); S12(4)=S12(4)+OME(i)*XJA*F2*(Cmat(6,1)*P11_1+Cmat(6,6)*P11_2+Cmat(6,6)*P 21_1+Cmat(6,2)*P21_2+Cmat(6,5)*P31_1+Cmat(6,4)*P31_2); S12(5)=S12(5)+OME(i)*XJA*F2*(Cmat(6,1)*P12_1+Cmat(6,6)*P12_2+Cmat(6,6)*P 22_1+Cmat(6,2)*P22_2+Cmat(6,5)*P32_1+Cmat(6,4)*P32_2); S12(6)=S12(6)+OME(i)*XJA*F3*(Cmat(6,1)*P13_1+Cmat(6,6)*P13_2+Cmat(6,6)*P 23_1+Cmat(6,2)*P23_2+Cmat(6,5)*P33_1+Cmat(6,4)*P33_2); S12(7)=S12(7)+OME(i)*XJA*F3*(Cmat(6,1)*P11_1+Cmat(6,6)*P11_2+Cmat(6,6)*P 21_1+Cmat(6,2)*P21_2+Cmat(6,5)*P31_1+Cmat(6,4)*P31_2); S12(8)=S12(8)+OME(i)*XJA*F3*(Cmat(6,1)*P12_1+Cmat(6,6)*P12_2+Cmat(6,6)*P 22_1+Cmat(6,2)*P22_2+Cmat(6,5)*P32_1+Cmat(6,4)*P32_2);
88 Chapter B. Codes S12(9)=S12(9)+OME(i)*XJA*F3*(Cmat(6,1)*P13_1+Cmat(6,6)*P13_2+Cmat(6,6)*P 23_1+Cmat(6,2)*P23_2+Cmat(6,5)*P33_1+Cmat(6,4)*P33_2); S13(1)=S13(1)+OME(i)*XJA*F1*(Cmat(5,1)*P11_1+Cmat(5,6)*P11_2+Cmat(5,6)*P 21_1+Cmat(5,2)*P21_2+Cmat(5,5)*P31_1+Cmat(5,4)*P31_2); S13(2)=S13(2)+OME(i)*XJA*F1*(Cmat(5,1)*P12_1+Cmat(5,6)*P12_2+Cmat(5,6)*P 22_1+Cmat(5,2)*P22_2+Cmat(5,5)*P32_1+Cmat(5,4)*P32_2); S13(3)=S13(3)+OME(i)*XJA*F1*(Cmat(5,1)*P13_1+Cmat(5,6)*P13_2+Cmat(5,6)*P 23_1+Cmat(5,2)*P23_2+Cmat(5,5)*P33_1+Cmat(5,4)*P33_2); S13(4)=S13(4)+OME(i)*XJA*F2*(Cmat(5,1)*P11_1+Cmat(5,6)*P11_2+Cmat(5,6)*P 21_1+Cmat(5,2)*P21_2+Cmat(5,5)*P31_1+Cmat(5,4)*P31_2); S13(5)=S13(5)+OME(i)*XJA*F2*(Cmat(5,1)*P12_1+Cmat(5,6)*P12_2+Cmat(5,6)*P 22_1+Cmat(5,2)*P22_2+Cmat(5,5)*P32_1+Cmat(5,4)*P32_2); S13(6)=S13(6)+OME(i)*XJA*F3*(Cmat(5,1)*P13_1+Cmat(5,6)*P13_2+Cmat(5,6)*P 23_1+Cmat(5,2)*P23_2+Cmat(5,5)*P33_1+Cmat(5,4)*P33_2); S13(7)=S13(7)+OME(i)*XJA*F3*(Cmat(5,1)*P11_1+Cmat(5,6)*P11_2+Cmat(5,6)*P 21_1+Cmat(5,2)*P21_2+Cmat(5,5)*P31_1+Cmat(5,4)*P31_2); S13(8)=S13(8)+OME(i)*XJA*F3*(Cmat(5,1)*P12_1+Cmat(5,6)*P12_2+Cmat(5,6)*P 22_1+Cmat(5,2)*P22_2+Cmat(5,5)*P32_1+Cmat(5,4)*P32_2); S13(9)=S13(9)+OME(i)*XJA*F3*(Cmat(5,1)*P13_1+Cmat(5,6)*P13_2+Cmat(5,6)*P 23_1+Cmat(5,2)*P23_2+Cmat(5,5)*P33_1+Cmat(5,4)*P33_2); S22(1)=S22(1)+OME(i)*XJA*F1*(Cmat(2,1)*P11_1+Cmat(2,6)*P11_2+Cmat(2,6)*P 21_1+Cmat(2,2)*P21_2+Cmat(2,5)*P31_1+Cmat(2,4)*P31_2); S22(2)=S22(2)+OME(i)*XJA*F1*(Cmat(2,1)*P12_1+Cmat(2,6)*P12_2+Cmat(2,6)*P 22_1+Cmat(2,2)*P22_2+Cmat(2,5)*P32_1+Cmat(2,4)*P32_2); S22(3)=S22(3)+OME(i)*XJA*F1*(Cmat(2,1)*P13_1+Cmat(2,6)*P13_2+Cmat(2,6)*P 23_1+Cmat(2,2)*P23_2+Cmat(2,5)*P33_1+Cmat(2,4)*P33_2); S22(4)=S22(4)+OME(i)*XJA*F2*(Cmat(2,1)*P11_1+Cmat(2,6)*P11_2+Cmat(2,6)*P 21_1+Cmat(2,2)*P21_2+Cmat(2,5)*P31_1+Cmat(2,4)*P31_2); S22(5)=S22(5)+OME(i)*XJA*F2*(Cmat(2,1)*P12_1+Cmat(2,6)*P12_2+Cmat(2,6)*P 22_1+Cmat(2,2)*P22_2+Cmat(2,5)*P32_1+Cmat(2,4)*P32_2); S22(6)=S22(6)+OME(i)*XJA*F2*(Cmat(2,1)*P13_1+Cmat(2,6)*P13_2+Cmat(2,6)*P 23_1+Cmat(2,2)*P23_2+Cmat(2,5)*P33_1+Cmat(2,4)*P33_2); S22(7)=S22(7)+OME(i)*XJA*F3*(Cmat(2,1)*P11_1+Cmat(2,6)*P11_2+Cmat(2,6)*P 21_1+Cmat(2,2)*P21_2+Cmat(2,5)*P31_1+Cmat(2,4)*P31_2); S22(8)=S22(8)+OME(i)*XJA*F3*(Cmat(2,1)*P12_1+Cmat(2,6)*P12_2+Cmat(2,6)*P 22_1+Cmat(2,2)*P22_2+Cmat(2,5)*P32_1+Cmat(2,4)*P32_2); S22(9)=S22(9)+OME(i)*XJA*F3*(Cmat(2,1)*P13_1+Cmat(2,6)*P13_2+Cmat(2,6)*P 23_1+Cmat(2,2)*P23_2+Cmat(2,5)*P33_1+Cmat(2,4)*P33_2); S23(1)=S23(1)+OME(i)*XJA*F1*(Cmat(4,1)*P11_1+Cmat(4,6)*P11_2+Cmat(4,6)*P 21_1+Cmat(4,2)*P21_2+Cmat(4,5)*P31_1+Cmat(4,4)*P31_2); S23(2)=S23(2)+OME(i)*XJA*F1*(Cmat(4,1)*P12_1+Cmat(4,6)*P12_2+Cmat(4,6)*P 22_1+Cmat(4,2)*P22_2+Cmat(4,5)*P32_1+Cmat(4,4)*P32_2); S23(3)=S23(3)+OME(i)*XJA*F1*(Cmat(4,1)*P13_1+Cmat(4,6)*P13_2+Cmat(4,6)*P 23_1+Cmat(4,2)*P23_2+Cmat(4,5)*P33_1+Cmat(4,4)*P33_2); S23(4)=S23(4)+OME(i)*XJA*F2*(Cmat(4,1)*P11_1+Cmat(4,6)*P11_2+Cmat(4,6)*P 21_1+Cmat(4,2)*P21_2+Cmat(4,5)*P31_1+Cmat(4,4)*P31_2); S23(5)=S23(5)+OME(i)*XJA*F2*(Cmat(4,1)*P12_1+Cmat(4,6)*P12_2+Cmat(4,6)*P 22_1+Cmat(4,2)*P22_2+Cmat(4,5)*P32_1+Cmat(4,4)*P32_2); S23(6)=S23(6)+OME(i)*XJA*F2*(Cmat(4,1)*P13_1+Cmat(4,6)*P13_2+Cmat(4,6)*P 23_1+Cmat(4,2)*P23_2+Cmat(4,5)*P33_1+Cmat(4,4)*P33_2); S23(7)=S23(7)+OME(i)*XJA*F3*(Cmat(4,1)*P11_1+Cmat(4,6)*P11_2+Cmat(4,6)*P 21_1+Cmat(4,2)*P21_2+Cmat(4,5)*P31_1+Cmat(4,4)*P31_2); S23(8)=S23(8)+OME(i)*XJA*F3*(Cmat(4,1)*P12_1+Cmat(4,6)*P12_2+Cmat(4,6)*P 22_1+Cmat(4,2)*P22_2+Cmat(4,5)*P32_1+Cmat(4,4)*P32_2);
B.9 SIGMAEQIGABEMANI 89 S23(9)=S23(9)+OME(i)*XJA*F3*(Cmat(4,1)*P13_1+Cmat(4,6)*P13_2+Cmat(4,6)*P 23_1+Cmat(4,2)*P23_2+Cmat(4,5)*P33_1+Cmat(4,4)*P33_2); S33(1)=S33(1)+OME(i)*XJA*F1*(Cmat(3,1)*P11_1+Cmat(3,6)*P11_2+Cmat(3,6)*P 21_1+Cmat(3,2)*P21_2+Cmat(3,5)*P31_1+Cmat(3,4)*P31_2); S33(2)=S33(2)+OME(i)*XJA*F1*(Cmat(3,1)*P12_1+Cmat(3,6)*P12_2+Cmat(3,6)*P 22_1+Cmat(3,2)*P22_2+Cmat(3,5)*P32_1+Cmat(3,4)*P32_2); S33(3)=S33(3)+OME(i)*XJA*F1*(Cmat(3,1)*P13_1+Cmat(3,6)*P13_2+Cmat(3,6)*P 23_1+Cmat(3,2)*P23_2+Cmat(3,5)*P33_1+Cmat(3,4)*P33_2); S33(4)=S33(4)+OME(i)*XJA*F2*(Cmat(3,1)*P11_1+Cmat(3,6)*P11_2+Cmat(3,6)*P 21_1+Cmat(3,2)*P21_2+Cmat(3,5)*P31_1+Cmat(3,4)*P31_2); S33(5)=S33(5)+OME(i)*XJA*F2*(Cmat(3,1)*P12_1+Cmat(3,6)*P12_2+Cmat(3,6)*P 22_1+Cmat(3,2)*P22_2+Cmat(3,5)*P32_1+Cmat(3,4)*P32_2); S33(6)=S33(6)+OME(i)*XJA*F2*(Cmat(3,1)*P13_1+Cmat(3,6)*P13_2+Cmat(3,6)*P 23_1+Cmat(3,2)*P23_2+Cmat(3,5)*P33_1+Cmat(3,4)*P33_2); S33(7)=S33(7)+OME(i)*XJA*F3*(Cmat(3,1)*P11_1+Cmat(3,6)*P11_2+Cmat(3,6)*P 21_1+Cmat(3,2)*P21_2+Cmat(3,5)*P31_1+Cmat(3,4)*P31_2); S33(8)=S33(8)+OME(i)*XJA*F3*(Cmat(3,1)*P12_1+Cmat(3,6)*P12_2+Cmat(3,6)*P 22_1+Cmat(3,2)*P22_2+Cmat(3,5)*P32_1+Cmat(3,4)*P32_2); S33(9)=S33(9)+OME(i)*XJA*F3*(Cmat(3,1)*P13_1+Cmat(3,6)*P13_2+Cmat(3,6)*P 23_1+Cmat(3,2)*P23_2+Cmat(3,5)*P33_1+Cmat(3,4)*P33_2); D11(1)=D11(1)+OME(i)*XJA*F1*(Cmat(1,1)*U11_1+Cmat(1,6)*U11_2+Cmat(1,6)*U 21_1+Cmat(1,2)*U21_2+Cmat(1,5)*U31_1+Cmat(1,4)*U31_2); D11(2)=D11(2)+OME(i)*XJA*F1*(Cmat(1,1)*U12_1+Cmat(1,6)*U12_2+Cmat(1,6)*U 22_1+Cmat(1,2)*U22_2+Cmat(1,5)*U32_1+Cmat(1,4)*U32_2); D11(3)=D11(3)+OME(i)*XJA*F1*(Cmat(1,1)*U13_1+Cmat(1,6)*U13_2+Cmat(1,6)*U 23_1+Cmat(1,2)*U23_2+Cmat(1,5)*U33_1+Cmat(1,4)*U33_2); D11(4)=D11(4)+OME(i)*XJA*F2*(Cmat(1,1)*U11_1+Cmat(1,6)*U11_2+Cmat(1,6)*U 21_1+Cmat(1,2)*U21_2+Cmat(1,5)*U31_1+Cmat(1,4)*U31_2); D11(5)=D11(5)+OME(i)*XJA*F2*(Cmat(1,1)*U12_1+Cmat(1,6)*U12_2+Cmat(1,6)*U 22_1+Cmat(1,2)*U22_2+Cmat(1,5)*U32_1+Cmat(1,4)*U32_2); D11(6)=D11(6)+OME(i)*XJA*F2*(Cmat(1,1)*U13_1+Cmat(1,6)*U13_2+Cmat(1,6)*U 23_1+Cmat(1,2)*U23_2+Cmat(1,5)*U33_1+Cmat(1,4)*U33_2); D11(7)=D11(7)+OME(i)*XJA*F3*(Cmat(1,1)*U11_1+Cmat(1,6)*U11_2+Cmat(1,6)*U 21_1+Cmat(1,2)*U21_2+Cmat(1,5)*U31_1+Cmat(1,4)*U31_2); D11(8)=D11(8)+OME(i)*XJA*F3*(Cmat(1,1)*U12_1+Cmat(1,6)*U12_2+Cmat(1,6)*U 22_1+Cmat(1,2)*U22_2+Cmat(1,5)*U32_1+Cmat(1,4)*U32_2); D11(9)=D11(9)+OME(i)*XJA*F3*(Cmat(1,1)*U13_1+Cmat(1,6)*U13_2+Cmat(1,6)*U 23_1+Cmat(1,2)*U23_2+Cmat(1,5)*U33_1+Cmat(1,4)*U33_2); D12(1)=D12(1)+OME(i)*XJA*F1*(Cmat(6,1)*U11_1+Cmat(6,6)*U11_2+Cmat(6,6)*U 21_1+Cmat(6,2)*U21_2+Cmat(6,5)*U31_1+Cmat(6,4)*U31_2); D12(2)=D12(2)+OME(i)*XJA*F1*(Cmat(6,1)*U12_1+Cmat(6,6)*U12_2+Cmat(6,6)*U 22_1+Cmat(6,2)*U22_2+Cmat(6,5)*U32_1+Cmat(6,4)*U32_2); D12(3)=D12(3)+OME(i)*XJA*F1*(Cmat(6,1)*U13_1+Cmat(6,6)*U13_2+Cmat(6,6)*U 23_1+Cmat(6,2)*U23_2+Cmat(6,5)*U33_1+Cmat(6,4)*U33_2); D12(4)=D12(4)+OME(i)*XJA*F2*(Cmat(6,1)*U11_1+Cmat(6,6)*U11_2+Cmat(6,6)*U 21_1+Cmat(6,2)*U21_2+Cmat(6,5)*U31_1+Cmat(6,4)*U31_2); D12(5)=D12(5)+OME(i)*XJA*F2*(Cmat(6,1)*U12_1+Cmat(6,6)*U12_2+Cmat(6,6)*U 22_1+Cmat(6,2)*U22_2+Cmat(6,5)*U32_1+Cmat(6,4)*U32_2); D12(6)=D12(6)+OME(i)*XJA*F2*(Cmat(6,1)*U13_1+Cmat(6,6)*U13_2+Cmat(6,6)*U 23_1+Cmat(6,2)*U23_2+Cmat(6,5)*U33_1+Cmat(6,4)*U33_2); D12(7)=D12(7)+OME(i)*XJA*F3*(Cmat(6,1)*U11_1+Cmat(6,6)*U11_2+Cmat(6,6)*U 21_1+Cmat(6,2)*U21_2+Cmat(6,5)*U31_1+Cmat(6,4)*U31_2); D12(8)=D12(8)+OME(i)*XJA*F3*(Cmat(6,1)*U12_1+Cmat(6,6)*U12_2+Cmat(6,6)*U 22_1+Cmat(6,2)*U22_2+Cmat(6,5)*U32_1+Cmat(6,4)*U32_2);
90 Chapter B. Codes D12(9)=D12(9)+OME(i)*XJA*F3*(Cmat(6,1)*U13_1+Cmat(6,6)*U13_2+Cmat(6,6)*U 23_1+Cmat(6,2)*U23_2+Cmat(6,5)*U33_1+Cmat(6,4)*U33_2); D13(1)=D13(1)+OME(i)*XJA*F1*(Cmat(5,1)*U11_1+Cmat(5,6)*U11_2+Cmat(5,6)*U 21_1+Cmat(5,2)*U21_2+Cmat(5,5)*U31_1+Cmat(5,4)*U31_2); D13(2)=D13(2)+OME(i)*XJA*F1*(Cmat(5,1)*U12_1+Cmat(5,6)*U12_2+Cmat(5,6)*U 22_1+Cmat(5,2)*U22_2+Cmat(5,5)*U32_1+Cmat(5,4)*U32_2); D13(3)=D13(3)+OME(i)*XJA*F1*(Cmat(5,1)*U13_1+Cmat(5,6)*U13_2+Cmat(5,6)*U 23_1+Cmat(5,2)*U23_2+Cmat(5,5)*U33_1+Cmat(5,4)*U33_2); D13(4)=D13(4)+OME(i)*XJA*F2*(Cmat(5,1)*U11_1+Cmat(5,6)*U11_2+Cmat(5,6)*U 21_1+Cmat(5,2)*U21_2+Cmat(5,5)*U31_1+Cmat(5,4)*U31_2); D13(5)=D13(5)+OME(i)*XJA*F2*(Cmat(5,1)*U12_1+Cmat(5,6)*U12_2+Cmat(5,6)*U 22_1+Cmat(5,2)*U22_2+Cmat(5,5)*U32_1+Cmat(5,4)*U32_2); D13(6)=D13(6)+OME(i)*XJA*F2*(Cmat(5,1)*U13_1+Cmat(5,6)*U13_2+Cmat(5,6)*U 23_1+Cmat(5,2)*U23_2+Cmat(5,5)*U33_1+Cmat(5,4)*U33_2); D13(7)=D13(7)+OME(i)*XJA*F3*(Cmat(5,1)*U11_1+Cmat(5,6)*U11_2+Cmat(5,6)*U 21_1+Cmat(5,2)*U21_2+Cmat(5,5)*U31_1+Cmat(5,4)*U31_2); D13(8)=D13(8)+OME(i)*XJA*F3*(Cmat(5,1)*U12_1+Cmat(5,6)*U12_2+Cmat(5,6)*U 22_1+Cmat(5,2)*U22_2+Cmat(5,5)*U32_1+Cmat(5,4)*U32_2); D13(9)=D13(9)+OME(i)*XJA*F3*(Cmat(5,1)*U13_1+Cmat(5,6)*U13_2+Cmat(5,6)*U 23_1+Cmat(5,2)*U23_2+Cmat(5,5)*U33_1+Cmat(5,4)*U33_2); D22(1)=D22(1)+OME(i)*XJA*F1*(Cmat(2,1)*U11_1+Cmat(2,6)*U11_2+Cmat(2,6)*U 21_1+Cmat(2,2)*U21_2+Cmat(2,5)*U31_1+Cmat(2,4)*U31_2); D22(2)=D22(2)+OME(i)*XJA*F1*(Cmat(2,1)*U12_1+Cmat(2,6)*U12_2+Cmat(2,6)*U 22_1+Cmat(2,2)*U22_2+Cmat(2,5)*U32_1+Cmat(2,4)*U32_2); D22(3)=D22(3)+OME(i)*XJA*F1*(Cmat(2,1)*U13_1+Cmat(2,6)*U13_2+Cmat(2,6)*U 23_1+Cmat(2,2)*U23_2+Cmat(2,5)*U33_1+Cmat(2,4)*U33_2); D22(4)=D22(4)+OME(i)*XJA*F2*(Cmat(2,1)*U11_1+Cmat(2,6)*U11_2+Cmat(2,6)*U 21_1+Cmat(2,2)*U21_2+Cmat(2,5)*U31_1+Cmat(2,4)*U31_2); D22(5)=D22(5)+OME(i)*XJA*F2*(Cmat(2,1)*U12_1+Cmat(2,6)*U12_2+Cmat(2,6)*U 22_1+Cmat(2,2)*U22_2+Cmat(2,5)*U32_1+Cmat(2,4)*U32_2); D22(6)=D22(6)+OME(i)*XJA*F2*(Cmat(2,1)*U13_1+Cmat(2,6)*U13_2+Cmat(2,6)*U 23_1+Cmat(2,2)*U23_2+Cmat(2,5)*U33_1+Cmat(2,4)*U33_2); D22(7)=D22(7)+OME(i)*XJA*F3*(Cmat(2,1)*U11_1+Cmat(2,6)*U11_2+Cmat(2,6)*U 21_1+Cmat(2,2)*U21_2+Cmat(2,5)*U31_1+Cmat(2,4)*U31_2); D22(8)=D22(8)+OME(i)*XJA*F3*(Cmat(2,1)*U12_1+Cmat(2,6)*U12_2+Cmat(2,6)*U 22_1+Cmat(2,2)*U22_2+Cmat(2,5)*U32_1+Cmat(2,4)*U32_2); D22(9)=D22(9)+OME(i)*XJA*F3*(Cmat(2,1)*U13_1+Cmat(2,6)*U13_2+Cmat(2,6)*U 23_1+Cmat(2,2)*U23_2+Cmat(2,5)*U33_1+Cmat(2,4)*U33_2); D23(1)=D23(1)+OME(i)*XJA*F1*(Cmat(4,1)*U11_1+Cmat(4,6)*U11_2+Cmat(4,6)*U 21_1+Cmat(4,2)*U21_2+Cmat(4,5)*U31_1+Cmat(4,4)*U31_2); D23(2)=D23(2)+OME(i)*XJA*F1*(Cmat(4,1)*U12_1+Cmat(4,6)*U12_2+Cmat(4,6)*U 22_1+Cmat(4,2)*U22_2+Cmat(4,5)*U32_1+Cmat(4,4)*U32_2); D23(3)=D23(3)+OME(i)*XJA*F1*(Cmat(4,1)*U13_1+Cmat(4,6)*U13_2+Cmat(4,6)*U 23_1+Cmat(4,2)*U23_2+Cmat(4,5)*U33_1+Cmat(4,4)*U33_2); D23(4)=D23(4)+OME(i)*XJA*F2*(Cmat(4,1)*U11_1+Cmat(4,6)*U11_2+Cmat(4,6)*U 21_1+Cmat(4,2)*U21_2+Cmat(4,5)*U31_1+Cmat(4,4)*U31_2); D23(5)=D23(5)+OME(i)*XJA*F2*(Cmat(4,1)*U12_1+Cmat(4,6)*U12_2+Cmat(4,6)*U 22_1+Cmat(4,2)*U22_2+Cmat(4,5)*U32_1+Cmat(4,4)*U32_2); D23(6)=D23(6)+OME(i)*XJA*F2*(Cmat(4,1)*U13_1+Cmat(4,6)*U13_2+Cmat(4,6)*U 23_1+Cmat(4,2)*U23_2+Cmat(4,5)*U33_1+Cmat(4,4)*U33_2); D23(7)=D23(7)+OME(i)*XJA*F3*(Cmat(4,1)*U11_1+Cmat(4,6)*U11_2+Cmat(4,6)*U 21_1+Cmat(4,2)*U21_2+Cmat(4,5)*U31_1+Cmat(4,4)*U31_2); D23(8)=D23(8)+OME(i)*XJA*F3*(Cmat(4,1)*U12_1+Cmat(4,6)*U12_2+Cmat(4,6)*U 22_1+Cmat(4,2)*U22_2+Cmat(4,5)*U32_1+Cmat(4,4)*U32_2);
B.9 SIGMAEQIGABEMANI 91 D23(9)=D23(9)+OME(i)*XJA*F3*(Cmat(4,1)*U13_1+Cmat(4,6)*U13_2+Cmat(4,6)*U 23_1+Cmat(4,2)*U23_2+Cmat(4,5)*U33_1+Cmat(4,4)*U33_2); D33(1)=D33(1)+OME(i)*XJA*F1*(Cmat(3,1)*U11_1+Cmat(3,6)*U11_2+Cmat(3,6)*U 21_1+Cmat(3,2)*U21_2+Cmat(3,5)*U31_1+Cmat(3,4)*U31_2); D33(2)=D33(2)+OME(i)*XJA*F1*(Cmat(3,1)*U12_1+Cmat(3,6)*U12_2+Cmat(3,6)*U 22_1+Cmat(3,2)*U22_2+Cmat(3,5)*U32_1+Cmat(3,4)*U32_2); D33(3)=D33(3)+OME(i)*XJA*F1*(Cmat(3,1)*U13_1+Cmat(3,6)*U13_2+Cmat(3,6)*U 23_1+Cmat(3,2)*U23_2+Cmat(3,5)*U33_1+Cmat(3,4)*U33_2); D33(4)=D33(4)+OME(i)*XJA*F2*(Cmat(3,1)*U11_1+Cmat(3,6)*U11_2+Cmat(3,6)*U 21_1+Cmat(3,2)*U21_2+Cmat(3,5)*U31_1+Cmat(3,4)*U31_2); D33(5)=D33(5)+OME(i)*XJA*F2*(Cmat(3,1)*U12_1+Cmat(3,6)*U12_2+Cmat(3,6)*U 22_1+Cmat(3,2)*U22_2+Cmat(3,5)*U32_1+Cmat(3,4)*U32_2); D33(6)=D33(6)+OME(i)*XJA*F2*(Cmat(3,1)*U13_1+Cmat(3,6)*U13_2+Cmat(3,6)*U 23_1+Cmat(3,2)*U23_2+Cmat(3,5)*U33_1+Cmat(3,4)*U33_2); D33(7)=D33(7)+OME(i)*XJA*F3*(Cmat(3,1)*U11_1+Cmat(3,6)*U11_2+Cmat(3,6)*U 21_1+Cmat(3,2)*U21_2+Cmat(3,5)*U31_1+Cmat(3,4)*U31_2); D33(8)=D33(8)+OME(i)*XJA*F3*(Cmat(3,1)*U12_1+Cmat(3,6)*U12_2+Cmat(3,6)*U 22_1+Cmat(3,2)*U22_2+Cmat(3,5)*U32_1+Cmat(3,4)*U32_2); D33(9)=D33(9)+OME(i)*XJA*F3*(Cmat(3,1)*U13_1+Cmat(3,6)*U13_2+Cmat(3,6)*U 23_1+Cmat(3,2)*U23_2+Cmat(3,5)*U33_1+Cmat(3,4)*U33_2); end D11=D11*JAC; D12=D12*JAC; D13=D13*JAC; D22=D22*JAC; D23=D23*JAC; D33=D33*JAC; S11=S11*JAC; S12=S12*JAC; S13=S13*JAC; S22=S22*JAC; S23=S23*JAC; S33=S33*JAC; end
References [1] B.A. Auld. Acoustic Fields and Waves in Solids. Number v. 1 in A Wiley-Interscience publication. Wiley, 1973. [2] Y. Bazilevs, V.M. Calo, J.A. Cottrell, J.A. Evans, T.J.R. Hughes, S. Lipton, M.A. Scott, and T.W. Sederberg. Isogeometric analysis using t-splines. Computer Methods in Applied Mechanics and Engineering, 199(5):229 – 263, 2010. Computational Geometry and Analysis. [3] Y. Bazilevs, J.A. Cottrell, and T.J.R. Hughes. Isogeometric Analysis. WILEY, Singapore, 2009. [4] C.A. Brebbia and J. Domínguez. Boundary elements. [accessed: 1st july 2015]. http://www.boundaryelements.com/index.php. [5] C.A. Brebbia and J. Dominguez. Boundary Elements An Introductory Course. WITPRESS, Southampton, 1996. [6] J.C. Marín A.Barroso F. París, J.Cañas. Introducción al análisis y diseño con materiales compuestos. Universidad de Sevilla. Escuela Técnica Superior de Ingenieros Industriales, 2008. [7] F. (2005) García Sánchez. Estudio numérico de problemas de fractura en materiales anisótropos elásticos y piezoeléctricos. [8] M. Guiggiani and P. Casalini. Direct computation of cauchy principal value integrals in advanced boundary elements. Journal for Numerical Methods in Engineering, 1987; 24: 1711-1720. [9] R.W. Johnson. Higher order b-spline collocation at the greville abcissae. Applied Numerical Mathematics, 2005; 52: 63-75. [10] Shivi Kesarwani. Polymer composites in aviation sector. International Journal of Engineering Research and Technology, V6, 06 2017. [11] Fred Kocks, Carlos Tomé, and H-R Wenk. Texture and Anisotropy. Preferred Orientations in Polycrystals and Their Effect on Material Properties. 01 2000. [12] Enayat Mahajerin and David L. Sikarskie. Boundary element study of a loaded hole in an orthotropic plate. Journal of Composite Materials, 20(4):375–389, 1986. [13] V. Mallardo, E. Roucco, J.Trevelyan, and T.Rabczuk. An improved isogeometric boundary element method approach in two dimensional elastostatic. CMES, 2014; 102: 373-391. [14] V. Mantič and F. París. Existence and evaluation of the two free terms in the hypersingular boundary integral equation of potential theory. Engineering Analysis with Boundary Elements, 16(3):253 – 260, 1995. [15] V. Mantič and F. París. Explicit formulae of the integral kernels and c-matrix in the 2d somigliana identity for orthotropic materials. Engineering Analysis with Boundary Elements, 15(3):283 – 288, 1995. 93
94 References [16] V. Mantič and F. París. Advanced formulation of the bies for elastic anisotropic materials in 2d. Advances in Boundary Element Techniques V, (3):57 – 62, 2004. [17] J.J. Pérez-Gavilán. Introducción a los elementos de frontera. 2006. [accessed: 1st july 2015]. [18] B. Simpson. Isogeometric-bem. 2012. [accessed: 14th september 2020]. https://github.com/bobbiesimpson/isogeometric-bem. [19] R.N. Simpson, J. Trevelyan, S.P.A. Bordas, and T. Rabczuk. A two-dimensional isogeometric boundary element method for elastostactic anaylisis. Computer Methods in Applied Mechanics and Engineering, 2012; 209-212: 87-100. [20] T.C.T.Ting. Anisotropic Elasticity Theory and Applications. 1996. [21] Ping Wang, Jinlan Xu, Jiansong Deng, and Falai Chen. Adaptive isogeometric analysis using rational pht-splines. Computer-Aided Design, 43(11):1438 – 1448, 2011. Solid and Physical Modeling 2011.