H_1hn2-wavelet Galerkin BEM and its application to the radiosity equation [Elektronische Ressource] / vorgelegt von Ulf Kähler

Description

2H -wavelet Galerkin BEM and itsapplication to the radiosity equationDISSERTATIONzur Erlangung des akademischen GradesDoctor rerum naturalium(Dr. rer. nat.)TECHNISCHE UNIVERSITAT CHEMNITZFakultat fur Mathematikvorgelegt von Dipl.-Math. Ulf Kahlergeb. am 28.07.1978 in StralsundBetreuer: Prof. Dr. Reinhold Schneider (TU Berlin)Gutachter: Prof. Dr. Dr. h.c. Wolfgang Hackbusch (MPI-MIS Leipzig)Prof. Dr. Arnd Meyer (TU Chemnitz)Prof. Dr. Reinhold Schneider (TU Berlin)Tag der Verteidigung: 5. November 2007Verfugb ar im MONARCH der TU Chemnitz:http://archiv.tu-chemnitz.de/pub/2007/0190ContentsIntroduction IIINotation index IX1 General de nitions and notations 11.1 Function spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Matrix notations and norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Relation notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 The boundary element method (BEM) 52.1 The geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 The boundary integral equation . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2.1 General considerations and assumptions . . . . . . . . . . . . . . . . . 62.2.2 The radiosity equation and similar problems. . . . . . . . . . . . . . . 72.2.3 The single layer potential operator . . . . . . . . . . . . . . . . . . . . 102.2.4 The double layer potential operator . . . . . . . . . . . . . . . .

Sujets

• Mathematics

Informations

• Publié par: technische_universitat_chemnitz
• Publié le: 01 janvier 2007
• Nombre de lectures: 8
• Langue: English
• Poids de l'ouvrage: 3 Mo

Extrait

2 Hwavelet Galerkin BEM application to the radiosity

D I S S E R T A T I O N

zur Erlangung des akademischen Grades Doctor rerum naturalium (Dr. rer. nat.)

and its equation

TECHNISCHEUNIVERSIT Ä TCHEMNITZ Fakultät für Mathematik

Betreuer: Gutachter:

vorgelegt von Dipl.Math. Ulf Kähler geb. am 28.07.1978 in Stralsund

Prof. Dr. Reinhold Schneider (TU Berlin) Prof. Dr. Dr. h.c. Wolfgang Hackbusch (MPIMIS Leipzig) Prof. Dr. Arnd Meyer (TU Chemnitz) Prof. Dr. Reinhold Schneider (TU Berlin)

Tag der Verteidigung: 5. November 2007

Verfügbar im MONARCH der TU Chemnitz: http://archiv.tuchemnitz.de/pub/2007/0190

Contents

Introduction

Notation index

1

2

3

4

General definitions and notations 1.1 Function spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Matrix notations and norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Relation notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

The boundary element method (BEM) 2.1 The geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The boundary integral equation . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 General considerations and assumptions . . . . . . . . . . . . . . . . . 2.2.2 The radiosity equation and similar problems . . . . . . . . . . . . . . . 2.2.3 The single layer potential operator . . . . . . . . . . . . . . . . . . . . 2.2.4 The double layer potential operator . . . . . . . . . . . . . . . . . . . 2.3 The Galerkin method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 The determination of the system matrix and the problem of quadratic costs . 2.5 Fast methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Fast multipole method . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Methods based on hierarchical matrices (H. . . . . . . . .matrices) . 2.5.3 Wavelet methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.4 The application of fast methods to the radiosity equation . . . . . . . 2.6 Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Multiscale basis – TauschWhite wavelets 3.1 The cluster tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Definition and classes of cluster trees . . . . . . . . . . . . . . . . . . . 3.1.2 Clustering methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Index notations for corresponding function sets . . . . . . . . . . . . . 3.2 Construction of the TauschWhite wavelets . . . . . . . . . . . . . . . . . . . 3.3 Wavelet configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Properties of the TauschWhite wavelets . . . . . . . . . . . . . . . . . . . . . 3.5 Enhancement of TauschWhite wavelets . . . . . . . . . . . . . . . . . . . . . 3.5.1 TauschWhite wavelets on binary trees . . . . . . . . . . . . . . . . . . 3.5.2 TauschWhite wavelets with vanishing moments with respect to the outer normals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Chapter conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Wavelet Galerkin BEM 4.1 Discrete wavelet transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Recursive determination of the system matrix . . . . . . . . . . . . . . . . . . 4.3 Apriori wavelet compression . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Assembly of the matrix pattern . . . . . . . . . . . . . . . . . . . . . . . . . .

I

III

IX

1 1 2 3

5 5 6 6 7 10 11 12 14 16 16 17 20 23 24

25 26 26 29 31 32 36 38 44 44

47 50

53 54 58 60 65

5

6

7

4.5 4.6 4.7

Aposteriori wavelet compression . . . . . . . . . . . . . . . . . . . . . . . . . Solving the linear system of equations . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 Introduction toHmatrices 5.1 Multidimensional interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . 2 5.2H. . . . . . . . . . . . . . . . . . . .approximation by interpolation . . . . 5.3 Hierarchical subdivision of the system matrix . . . . . . . . . . . . . . . . . . 5.4 Uniform matrix and matrixvector multiplication . . . . . . . . . . . . . . . . 5.5 Nested cluster basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 General estimates and conclusion . . . . . . . . . . . . . . . . . . . . . . . . .

2 HWavelet Galerkin method 6.1 Approximation of a single entry . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Nested scaling and wavelet cluster basis . . . . . . . . . . . . . . . . . . . . . 6.3 Fast recursive determination of matrix blocks . . . . . . . . . . . . . . . . . . 6.4 Fast algorithm for the computation of the system matrix . . . . . . . . . . . . 6.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Wavelet Radiosity 2 H. . . . . . . . . . . . . .wavelet Galerkin method for the unoccluded case Problems and geometrical considerations of the occluded case . . . . . . . . . 7.2.1 General geometrical considerations and assumptions . . . . . . . . . . 7.2.2 Visibility between panels . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.3 Visibility between clusters . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.4 The effect of partial visibility on fast methods . . . . . . . . . . . . . . Wavelet compression for the radiosity equation . . . . . . . . . . . . . . . . . 7.3.1 General approach to the wavelet compression for the different visibility cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.2 Estimates for the partial visibility case A . . . . . . . . . . . . . . . . 7.3.3 Estimates for the partial visibility case B . . . . . . . . . . . . . . . . 7.3.4 Estimates for the partial visibility case C . . . . . . . . . . . . . . . . 7.3.5 Estimates for the partial visibility case D . . . . . . . . . . . . . . . . 7.3.6 Final wavelet compression for the radiosity equation . . . . . . . . . . 2 Adaption of theHwavelet method onto the radiosity equation . . . . . . . . Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7.1 7.2 7.3 7.4 7.5 7.6

Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Bibliography

Theses

II

67 68 68

71 72 73 76 79 81 84

87 88 92 96 102 114 120

121 122 127 127 130 133 135 136

136 139 146 149 153 153 155 158 164

167

171

Introduction

In modern science and engineering a large amount of problems can be formulated as boundary integral equations. These equations can arise, e.g. from a reformulation of a partial differential equation in an exterior boundary value problem. In this thesis we consider among others the radiosity equation which is directly derived from a physical state. It is an important approach for the simulation of illumination in virtual scenes within computer graphics [3, 20]. Its solution provides a measure for the brightness of an object and takes the smoothness of shadows into account. Furthermore, the numerical treatment of the radiosity equation is similar to the treatment of the heat radiation equation. This means that all approaches which we will make to the radiosity equation in this thesis can also be applied to the heat radiation equation.

In general, boundary integral equations, i.e. Z λρ(x)−k(x, y)ρ(y)dΓ y =f(x), x∈Γ, Γ are solved numerically by the boundary element method (BEM). It uses a finite discreti sation of the boundary Γ similar to the finite element method (FEM). However, contrary 2 to the FEM a traditional discretisation leads to an associated system matrix withO(N) entries, whereNdescribes the number of degrees of freedom as well as the number of equations. There are several methods to reduce the resulting unacceptable quadratic costs in computation time and memory usage to linearlogarithmic or even linear costs like the multipole method[38, 54] orpanel clustering[43, 58]. The research on the last method has led to the important class of fast methods based onhierarchical matrices, which includes 2 Hmatrices[40, 41],Hmatrices[42, 35, 9],ACA[4] andHCA[10].

An alternative approach are the wavelet methods introduced in [5] and improved in [26, 27, 28, 29, 44, 56, 61]. For that, the Galerkin discretisation within the BEM is performed by means of a multiscale wavelet basis with a sufficient amount of vanishing moments, i.e., the wavelets are orthogonal to the traces of polynomials up to a certain degree. As a consequence, the resulting system matrix is quasisparse. This means that the matrix in the wavelet basis has onlyO(Nlog(N)) relevant entries for general geometries andO(N) relevant entries for smooth geometries. By use of a cutoff parameter [61] the matrix can apriori be compressed to a sparse system matrix, which provides the desired linearlogarithmic costs. Moreover, with a second aposteriori compression [44] the amount of relevant entries can be further reduced.

Unfortunately, special boundary integral operators cause special difficulties. Many integral operators possess singularities in their kernel functions, which are well studied. The radiosity equation additionally contains the socalled visibility function which causes further discontinuities. There are several approaches to handle these difficulties e.g. multipole methods [1, 2], wavelet approaches [20, 19, 63] or quasi montecarlo methods [49, 50]. Since the discontinuities destroy the structure on which the efficiency of hierarchical matrices and multipole methods are based, they prevent these methods from having linearlogarithmic costs for general geometries. The wavelet approach, calledwavelet radiosity, was proclaimed

III

to keep its linearlogarithmic costs. However, the discontinuities also prevent the usage of the standard cutoff parameter. The necessary new cutoff parameter for the special situation of the radiosity equation as well as the corresponding proofs have been missing in the literature until now. In this thesis we will close this gap for the 3dimensional radiosity equation. We will present the required new cutoff parameter and the corresponding proofs for a wide class 2 of geometries, which will allow us to compress the system matrix toO(NlogNFor) entries. that, we will use properties of wavelets concerning their multiscale structure presented in [61].

Independently whether we consider the radiosity equation or a general boundary integral equation, a fast implementation of the wavelet method is necessary for any application. However, it is not trivial to construct a wavelet method whose complexity, in terms of com putation and memory requirements, is proportional to the number of relevant entries in the system matrix. The wavelet method consists of four different parts: 1. the determination of the righthand side vector,

2. the computation of the system matrix,

3. the solution of the resulting linear system of equations and

4. the transfer of the solution vector into a singlescale basis for further applications. The first as well as the last part are well studied in form of the discrete wavelet transform. The third part is supported by the sparse system matrix and the presence of the wavelets. It remains to compute the system matrix in wavelet coordinates. This is the most difficult task of the wavelet method. Due to the multiscale structure of a wavelet its support can cover large parts of the boundary. Although onlyO(N) orO(Nlog(N)) entries have to be α computed, the cost for the computation of every single entry has to beO(1) orO(log (N)), even if the corresponding integration area covers nearly the whole boundary.

In the literature there are only a few effective wavelet methods which satisfy these demands. In [44] such an effective method has been presented for geometries which are sufficiently smooth and are given by a parameterisation. It computes the sparse system matrix with O(N) entries inO(N) operations using biorthogonal wavelets. However, for the present thesis we want to assume that the geometry is only given by its polygonal approximation. It is the representation of the geometry with the lowest amount of information, since we have only a set of triangles and the set of the corresponding points. As a consequence, we can only guarantee a compression toO(Nlog(N)) entries for the resulting system matrix due to the missing smoothness conditions. Furthermore, we can neither use the biorthogonal wavelets nor the method in [44] for the computation of the system matrix.

Therefore, we will use a multiscale basis which is called TauschWhite wavelet basis in the literature [65]. Its basis functions are recursively defined. They are a linear combination of the well known piecewise constant ansatz functions and provide the required vanishing moment condition. The necessary properties concerning the wavelet compression are proved in [45] following [65, 61]. Moreover, the recursive representation of the TauschWhite wavelets carries over to the entries of the system matrix. In order to solve the problem of the large supports of the wavelets, we use the ideas of another fast method. Among others 2 theHmatrix provides lowrank approximations of submatrices of the system matrix in a singlescale basis and parts of their hierarchical structure perfectly fit to the structure of 2 the TauschWhite wavelets. Therefore, we will combine parts of theHmatrices with the 2 recursively defined wavelets and construct a new wavelet method, which we callHwavelet Galerkin method.

2 This new method will requireO(Nlog(N)k) operations and memory efforts of O(N(log(N) +k)) to set up the sparse system matrix in the wavelet basis withO(Nlog(N))

IV

β 2 entries, wherek∼logNis the lowrank of the approximation used by theHmatrices. Once the system matrix is established, we only need furtherO(Nlog(N)) operations to solve the wavelet Galerkin scheme for any righthand side. This provides a great advantage for solutions for a huge amount of different righthand sides.

As mentioned before, the radiosity equation causes a special difficulty. The equation is written as Z hnˇ x , x−yihˇn y , y−xi B(x) =E(x) +r d (x)B(y)V is(x, y)dΓ y , x∈Γ =∂Ω, 4 πkx−yk Γ with the inner normalsnˇ x andnˇ y , the radiosityB(x), which describes the brightness of the surface, the emitted radiationE(xa light source, and the diffuse reflection coefficient), i.e. r d (xA short) describing the amount of radiation which is not absorbed by the surface. derivation of this equation is presented in subsection 2.2.2.

The problems concerning the radiosity equation arise from the visibility functionV is(x, y). It is defined by  1 if (x+t(y−x))∈allΩ for t∈(0,1), V is(x, y) := 0 elsewhere. The first difficulty of the visibility function is its determination for different surface parts. 3 2 A naive test for a discretisation with panelsNwould requireO(NFor) operations. N couples of panels we would have to test whether the remaining (N−2) panels blockade 2 their line of view. Depending on the geometry this costs can be reduced toO(N) or even O(Nlog(Nin this thesis we are more interested in a second problem. Therefore,)). However, we assume to have a suitable visibility test.

x µ y

V ′ y P ′ x I

Figure 1: The effect of visibility (left) and the system in a singlescale basis on a Lshape. (Blue represents nonzero entries).

matrix (right) of the radiosity equation zero entries, while nonblue represents

The second difficulty concerning the visibility function consists in its jumps or discontinuities, respectively. For a single pointxthe visibility function divides the geometry in a visible and a nonvisible part, cf. Figure 1 (left). For a whole partµof the surface situated between the pointsxandythe geometry is divided into a visible partV, a nonvisible partIand a part Pwhich is partly visible byµand contains the jumps from visible to nonvisible. In the Figure 1 (right) we can see the effect of the visibility function on the system matrix. Blocks which are completely blue correspond to nonvisible parts. Blocks which are yellow and red correspond to visible parts. Last but not least, blocks which are in parts blue and in parts red and yellow correspond to partly visible parts and contain the discontinuities.

V

Fast methods like the methods based on hierarchical matrices or the multipole method approximate blocks of the system matrix. For that, they require a certain smoothness of the kernel of the integral equation. As a consequence, the approximation can only be performed for visible and nonvisible blocks. All blocks which correspond to partly visible surface parts cannot be approximated. Hence, the complexity of time and memory usage cannot be reduced to linearlogarithmic. However, it is still possible for certain geometries to reduce 3 the costs toO(N), cf. [2]. 2

As mentioned above the wavelet methods suffer the same problem, since the wavelet compression is based on series expansions. The wavelet compression from [61] can only be applied to couples of wavelets whose supports are visible or nonvisible. With the new cutoff parameter, which we will present in this thesis, we can also provide a wavelet compression for wavelets whose supports are partly visible.

2 With a slight modification we will apply theHwavelet Galerkin method to the radiosity 2 equation. This method is based on theHmatrices. Therefore, we have time cost for the 3 2 β assembly of the system matrix in wavelet coordinates ofO(N k) withk∼logNsimilar 2 2 to theHdue to the new cutoff parameter we will set up a systemmatrices. However, 2 matrix with onlyO(NlogN) entries. With it, we can reduce the temporary memory cost 2 toO(N(logN+k)).

2 Similar to theHwavelet Galerkin method for general boundary integral equations, we 2 only need furtherO(NlogN) operations to solve the wavelet Galerkin scheme for any righthand side. This provides a great advantage as soon as the system matrix is assembled. Since the the visibility function’s information will be assimilated by the system matrix in wavelet coordinates, the visibility function has no further effect to the cost of the solution of further righthand sides. This is an even greater advantage since the righthand side of the radiosity equation represents the different light sources, while the system matrix contains the information of the geometry. E.g., the system matrix has to be set up once for the computation of the illumination of a room. Then, one can test the effect of various differ 2 ent strengths and positions of several lamps with time and memory costs of onlyO(NlogN).

We will tackle the mentioned tasks as follows. In the first short chapter we introduce some necessary notations and definitions.

The second chapter provides a compact introduction to the boundary element method, which includes a summary of the most important fast methods. In this context we formulate the general assumptions to the boundary integral equations as well as to the given geometry. Furthermore, we give a short derivation for the radiosity equation and finish the chapter with the formulation of the major goals of this thesis.

The third chapter is devoted to the TauschWhite wavelet basis. We will introduce the necessary hierarchical subdivision of the geometry called cluster tree and construct the TauschWhite wavelets on this tree. Following [45], we will present important estimates and properties of the TauschWhite wavelets. Besides, we will perform certain modifications to the wavelet functions which are necessary for further considerations.

The fourth chapter will deal with properties and estimates for a general wavelet Galerkin method. This includes the discrete wavelet transform, the presentation of the apriori compression as well as the aposteriori compression and remarks for the solution of the resulting linear system of equations. It also contains the first steps for the recursive computation of the matrix entries as well as an algorithm for the apriori determination of the matrix pattern of the compressed system matrix. Both will be used in chapter 6 by the

VI

2 fastHwavelet Galerkin method.

2 The fifth chapter provides an introduction toHwill shortly present thematrices. We 2 method and make necessary considerations for the upcoming construction of theHwavelet 2 Galerkin method which includes theHapproximation by interpolation and the nested cluster basis.

In the sixth chapter we will harvest the fruits of the previous chapters and will complete 2 theHwavelet Galerkin method with the algorithm for the fast computation of the system 2 matrix. For that, we will adapt theHapproximations for submatrices concerning the piecewise constant ansatz functions to the entries concerning the wavelet basis. Throughout this chapter we will combine this approach with the recursive computation of the entries, which provides an algorithm for the fast assembly of the system matrix. Numerical results are presented that support the complexity and error estimates.

The seventh and last chapter is devoted to the fast solution of the radiosity equation. In this chapter we will develop a new cutoff parameter for the wavelet compression for the radiosity equation concerning the cases in which the visibility function becomes effective and 2 apply a modifiedHFor that, we willwavelet Galerkin method to the radiosity equation. make certain assumptions to the geometry and specify the cases of jumps within the visibility functions. For these cases we will prove the necessary error estimates and estimates for the 2 amount of relevant entries. For the application of theHwavelet Galerkin method we adapt 2 theHFinally, we present numericalapproximations to the needs of the radiosity equation. results that support the mentioned complexity and error estimates.

VII

VIII

Notation

Function spaces

k k k C(Ω), C(Γ), C(Ω), pw k C(Γ) pw k,µ k,µ C(Ω), C(Γ), k,µ k,µ C pw (Ω), C pw (Γ) k,µ , k,µ k C C , C , C pw pw

q q q H , H(Γ), H(Ω) 2 2 2 L , L(Γ), L(Ω) ∞ ∞ ∞ L , L(Γ), L(Ω)

index

spaces of ktimes continuously or piecewise continuously differen tiable functions spaces of ktimes Hölder continuously or ktimes piecewise Hölder continuously diff. functions set of all surfaces with the corresponding smoothness or piecewise smoothness q q Sobolev spaces (H=H(Γ)) 2 2 spaces of quadratic Lebesque integrable functions (L=L(Γ)) ∞ ∞ spaces of functions with finite essential supremum (L=L(Γ))

Norms and scalar products

k ∙ k,k ∙ k 2 2 L(Ω)L(Γ) h∙,∙i,h∙,∙i 2 2 L(Ω)L(Γ) k ∙ k L(Ω) ,k ∙ k L(Γ) ∞ ∞ k ∙ k H s k ∙ k 1 ,k ∙ k 2 ,k ∙ k ∞ h∙,∙i

2 Lnorm in Ω or on Γ, resp., (k ∙ k 2 =k ∙ k 2 ) L L(Γ) inner product in Ω or on Γ, resp., (h∙,∙i 2 =h∙,∙i 2 ) L L(Γ) essential supremum in Ω or on Γ, resp., (k ∙ k L =k ∙ k L(Γ) ) ∞ ∞ s associated norm to the Sobolev spaceH(Γ) 1, spectral and infinity matrix norm n+1 usual scalar product inR

Operators and kernel functions

V,K k S (x, y), k D (x, y) R k R (x, y) V is(x, y) k Rad (x, y)

single or double layer potential operator kernel of the single or double layer potential operator operator of the radiosity equation kernel function of the radiosity equation visibility function of the radiosity equation kernel of the radiosity equation,k Rad (x, y) =k R (x, y)V is(x, y)

Geometrical objects

Γ N N h π i n i ,ˇn i µ, ν son pa ν , ν , νˆ i j ν J T I N ,I ν

polygonal approximation of the boundary number of panels step width simple panel outer or inner normal vector of the panelπ i clusters son, father or root cluster, resp. level concerning the clusterν finest level of the clustertree index set of all panels of Γ N or from the clusterν, resp.

IX