BOUNDARY ELEMENT METHOD MODYFICATIONS FOR USE IN SOME IMPEDANCE AND OPTICAL TOMOGRAPHY APPLICATIONS

The article presents two elements associated with the practice of application of the boundary element method. The first is associated with BEM ability to analyze an open boundary objects and application of infinite boundary elements in the area of mammography. The second element is associated with the damped wall investigations. Wall humidity and moisture represents heterogeneous environment (Functionally Graded Materials) which has to be treated in a special way.


Introduction
One of the advantages of the Boundary Element Method is its ability to handle open boundary problems. It is useful either when it is necessary to analyze outer field distribution or when there are some problems describing object surface and related boundary conditions. Another important tomograpic task is to find internal properties spatial distribution within non-homogenous region.
A good example of open boundary and related infinite boundary element is an application for early detection or screening examination of breast cancer using Optical Tomography. Using infinite boundary elements to create an open boundary breast model results in the significant reduction of the mesh size and allows us to avoid problems with unknown boundary conditions. In consequence, hardware requirements are much lower and the results are received faster without accuracy losses. The implementation of two main lines of infinite boundary elements will be discussed on the example of a typical hemispherical breast model.
Boundary Element Method (BEM) is employed in breast cancer investigations to find a forward problem solution and in turn to receive the internal image of breast tissue solving the reverse problem, as it is not possible to place detectors or to precisely define boundary conditions on the surface between breast and chest.
The traditional method used to avoid errors near the breast and chest boundary is to extend the mesh outside the region of interest, deep into the chest and to truncate it at some distance from the investigated breast tissue. Wrong boundary conditions or improper placement of such artificial boundary can introduce an unknown error if the truncation occurs too near. On the other hand, excessive mesh received while setting that boundary too far, increases the number of boundary elements and decreases the computational efficiency. One of the alternative solutions is to use infinite boundary elements and to create an open boundary model. The implementation of these elements can reduce the mesh and help avoid the problem of incorrect boundary conditions assumption.
Moreover, despite constant hardware development the process of finding the inner breast image in computer tomography is still a time-consuming task. The implementation of infinite boundary elements is an attractive option to solve these issues. Results can be found faster due to the significant reduction of the mesh size and without accuracy losses.
There are two main lines of infinite boundary elements development: decay functions infinite elements and mapped infinite elements. The first type uses special decay functions in conjunction with ordinary boundary element interpolation functions [2,3,5,23]. In that case the field variable tends monotonically to its far field value while reaching the element boundary adjacent to the infinite surroundings. Consequently, along finite element length, the variable changes in the way that reflects the physics of the problem up to infinity. The second type transforms the element from finite to infinite domain. The field variable will reach its far value following geometrical coordinates which extend into infinity [5,12,25]. Changing the basic interpolation functions interferes in the Boundary Element Method fundamental rules. The application of infinite elements in the open edge of the object usually requires using special quadratic boundary elements. This constitutes another complication for most mesh generators. On the other hand, the process of infinite elements incorporation into BEM is quite logical and results in significant hardware requirements and calculation time reduction. All major aspects of the implementation of both infinite elements types into BEM will be presented on the example of optical mammography.
Second presented subjectapplication of Functionally Graded Materials into BEM do not allow us to decrease BEM mesh size, like infinite elements mentioned above. Substituting multilayer BEM model related to non-homogenous environment by single layer model with applied FGM keeps full advantage of an important BEM feature, ie. reduction of model size.
Serious problem form many brick masonry historic buildings are their excessive moisture. Too much moisture causes a reduction in the compressive strength of both brick and mortar and to reduce the durability of the walls. Water, usually containing aggressive chemicals (chlorides, sulfates and nitrates), is transmitted to higher altitudes of the wall. High humidity of the walls encourages growth of fungi and mold, which have a negative impact on the human health. Humidity distribution, the type and concentration of salt in the wall are helpful in assessing the causes of damage and should be the basis for selecting the appropriate security methods related to excessive humidity and salinity.
Typical method used to find humidity distribution performed by drilling, sampling and checking samples moisture cannot be applied for continuous monitoring of dehumidification process. Relatively inexpensive method based on Electrical Impedance Tomography (EIT) [4,6,7,9,15,18] can be implemented. The wall conductivity is changing with humidity of the wall and its salinity, so object propertiesconductivity -varies smoothly with one ore more co-ordinates. Boundary Element Method (BEM) [10,17] was used to solve the EIT forward problem. BEM implementation requites to use either multiregional model or Functionally Graded Materials (FGM) to find spatial distribution of wall humidity.
Some theoretical aspects of BEM and FGM usage in wall dampness investigations [4], are within the scope of this paper.

Infinite boundary elements application
for mammography

Theory of infinite elements
The basic idea of the decay function infinite boundary element construction is that the standard basis interpolation functions Ni are multiplied by so called decay functions Di [2,3,5,23]. Two types of decay functions will be considered: reciprocal and exponential. Reciprocal decay functions for the decay in positive direction of ξ are as follows [5]: where i corresponds to the node number, (ξo ,ηo) is some origin point. This point must be outside the infinite element on the opposite side to the one which extends to infinity, n has to be greater than the highest power of ξ encountered in Ni . If the decay is in the positive ξ direction then ξo < −1. Respectively ηo < −1 for the decay in the positive η direction. This avoids a singularity within the element. For second order, eight-node quadrilateral elements, basis interpolation functions n=3 was chosen. Decay function infinite basis interpolation function becomes: For the exponential decay function of the form: where L[m] is the length which determines the severity of the decay, basis infinite interpolation functions are given by: The basis interpolation function Ni for standard eight-node quadrilateral isoparametric boundary elements are given by the following formulas (5): Decay function infinite elements based on 8 nodes, the second type, quadrilateral standard boundary elements are presented in figure 1. It should be noticed that it consists of 8 nodes. The only important thing is to keep the correct relation between decay function and node numbers, which decides in which direction the element, or more exactly the element properties (geometry remains unchanged) extends into infinity.
Mapped infinite boundary elements based on 8 nodes, the second type, quadrilateral standard boundary elements are presented in figure 2. It should be noticed that in this case the element will consist only of five nodes. Nodes 2, 3, 4 tend to infinity and will not take part in the calculations.
Despite its name the procedure for deriving these functions is quite logical and clearly described by Bettess [5, in chapter 4]. The infinite basis interpolation functions i M grow without limit as a coordinate approaches infinity, and are applied to the geometry. The ordinary basis interpolation functions Ni are applied to the field variables [5,12,23].
It is necessary to use these infinite basis interpolation functions to calculate the Jacobian and regularization transformations. The Jacobian is related to the transformation from two-dimensional intrinsic mesh structures to the global three dimensional coordinate system. In other words, we need to define the way in which we can pass from the x,y,z global Cartesian system to the ξ,η,ζ system defined over the element, where ξ and η are oblique coordinates and ζ is in the direction of the normal. The transformation for a given function Φ is related through the following: where the square matrix is the Jacobian matrix (or Jacoby matrix). IAPGOŚ 1/2017 p-ISSN 2083-0157, e-ISSN 2391-6761 A transformation of this type allows us to describe differentials of surface in the Cartesian system in terms of the curvilinear coordinates. A differential of area will be given by: and n is a normal vector to the surface element.
The coordinates x, y, z of a point on an element with the intrinsic coordinates ξ and η are given by: where: M i is a basis interpolation function at node i, x i , y i , z i represent coordinates of i-th node of the element. The required components of the unit outward normal (n x , n y , n z ) from Equation (8) are given by: The analyzed area as well as the corresponding integral equation will consist of both parts: finite and infinite surrounding of an open edge.
For debugging purposes, in case of ordinary basis interpolation functions one has to check if the sum of all basis interpolation functions is unity and if the sum of all their derivatives is zero. A simple test checks if each function has unit value at its own node and zero at the other nodes. For decay type functions and for nodes remaining in the calculations in case of mapping functions that condition is also fulfilled. There is no exact analogy for the nodes which escape into infinity at mapped infinity elements. Further tests using Zienkiewicz type of mapped infinite elements [25] are devised by Bettess [5].

Breast models descriptions
Three simple theoretical models of human breast were investigated. For all models one placement of the light source was presented -located near the bottom of the hemisphere model. The first model presented in fig.3 corresponds to the hemisphere with an additional cylindrical part at the bottom. The next two open boundary models use infinite boundary elements instead of the additional cylindrical part. The area of interest is limited to the hemisphere. The extra cylindrical part or the infinite elements rings are necessary to be added just to avoid possible errors, which would occur in case of mesh truncation at the bottom of the hemisphere. The standard boundary element model (Fig. 3) was constructed from 1536 second order, eightnode quadrilateral boundary elements and 4610 nodes. Half of the elements cover the hemisphere. The open boundary model consists of 832 standard boundary elements and 64 infinite elements based on eight-node second order quadrilateral boundary elements [12]. The number of nodes was reduced to 2625 nodes in the model with incorporated mapped infinite boundary elements and to 2753 nodes in the decay function infinite element usage.
The governing equation for the problem is the diffusion approximation of the transport equation [1,16,21,24] (Helmholtz (11) -assuming scattering and absorption are homogeneous). There are Robin boundary conditions (12) on surfaces [1,16,24]. In Diffusive Optical Tomography the distribution of the absorbing coefficient a  and the reduced scattering coefficient Where Φ stands for photon density, with different coefficients for breast tissue and for skeletal muscles on the basis [1,3,25] imposed. In the analyzed example the following breast tissue properties were taken [1,2,16]: The relevant boundary integral equation for surfaces covered by standard and infinite elements can be written as: n is the number of these sources, Φ stands for the photon density and G is the fundamental solution for the diffusion equation [1,16,21,24]. In 3D space for the diffusion equation the fundamental solution is [16]: The normal derivative of the Green function in the direction n can be written as:

Results and discussion
The values of /n were already presented above in figures 3, 4 and 5. It is difficult to compare these color maps, especially as each model and its solution have its own range of values, so the same color on all maps does not correspond to the same value. Taking that into account standard graphs were used to compare the results precisely. The values of /n module and the phase of the light at nodes lying on hemisphere circumference cross-section for y=0 are presented in figures 7 and 9 respectively. The values in nodes are presented in relation to angle Ψ (Fig. 6). To estimate the solution differences, the model with the extended bottom part (Fig. 3) was compared to those with infinite boundary elements implemented ( Fig. 4 and 5).
As it was mentioned in the introduction the implementation of infinite boundary elements requires changing one of its fundamental elementsthe basic interpolation functions. So, to present and to compare the different types of infinite boundary elements it was necessary to write our own BEM software in C++. The software uses some mathematical libraries from TOAST package [22]. Generally, all results are almost identical. Differences for module and phase are presented in figures 8 and 10 respectively. Only the part far from the light source shows some fluctuations, less than 1%, in results but it should be noticed that in that region the module is 10 orders of magnitude lower than near the source.

Laboratory model and problem formulation
Tested wall parameters, measurements and theoretical aspects of EIT will be presented. Real, prepared for the experiment brick wall presented in figure 11 was 1m height, 1m long and 0.51m thick. It was flooded up to half height by water. The measurements were made a day after the water was drained.
Coefficients a00, a10, a01, a11 and a2 where calculated in inverse problem solution in accordance with data measured in Impedance Tomography.
In a simple BEM usage fundamental solutions are known and tabulated. For Functionally Graded Materials it is necessary to find the adequate Green's function [8,11,13,14,19,20].
Green's function in two dimensions is derived for graded material represented by damped wall. The corresponding boundary integral equation formulations for this problem is derived, and the two-dimensional case is solved numerically.

Theoretical considerations
Boundary Integral Equation (BIE) [10,17] corresponding to the problem is: Comparing the above equation to the similar BIE for homogenous materials [10,17] two new elements can be found. At first (r) in equation (2) represents conductivity spatial distribution and secondly G* stands for so called modified Green function for non-homogenous materials. Curve Γ corresponds to the wall boundaries.
Modified Green function G* is not a Laplace equation fundamental solution. The differential equation for a potential function  defined on a surface S bounded by a wall boundary curve Γ, with an outward normal n, can be written as: The kernels are based on the Green's function G* defined as: where: r and r' are load and source points respectively. Following E. Kurgan [11] or A. Sutradhar and G. H. Paulino [13,14,19,20] and assuming that conductivity will fulfil a condition: we will receive applicable solution: where G is a typical Green's function for Laplace equation [10,17]. Matrices elements A and B [10,17] can be written as follows:

Mammography
All three types of infinite boundary elements: mapped, reciprocal decay functions and exponential decay functions offer almost identical results, similar to these achieved by using the standard model, which contains only finite boundary elements. The role of the infinite elements is to receive the correct solution in the area of interest, so all additional parts like the bottom cylinder in the standard model (Fig. 3) as well as the ring built from infinite elements around the hemisphere (Fig. 4 and 5) are neglected.
The advantages of using infinite elements are to avoid incorrect and unknown boundary conditions such as those on the surface between breast and chest in the breast models, to shorten the calculation time and to keep the accuracy similar to the standard solution (based on mesh extension, in our case mesh outside the hemisphere). All extended models were built from the same number of standard 8-node, second order, quadrilateral boundary elements. Mesh density on the additional surface related to the cylindrical part of the model was lower than on the hemisphere surface. This is a typical practical solution, as the additional part represents the region outside the zone of interests.
Reducing the number of mesh elements almost to 50% is fundamental to the inverse problem solution when the forward problem has to be calculated many times. The implementation of infinite boundary elements into the boundary element method improves the computational efficiency in the breast model almost 4 times and enables us to avoid the problem of setting incorrect boundary conditions.
The forward problem calculations took 4 minutes and 47 seconds in the case of 4160-node model. The model with infinite elements built, in total, from 2625 nodes, with incorporated mapped infinite boundary elements, and the model built from 2753 nodes, with the decay function infinite elements, took respectively 1 minute and 24 seconds and 1 minute and 30 seconds.. The reverse problem solution used to find the inner image of the analyzed objects requires calculating the forward problem many times.
The process of incorporating infinite elements into BEM calculation scheme is quite logical and generally related to the incorporation of new infinite basic approximation function and all further steps for its implications.
The main disadvantage of using infinite boundary elements is that not all mesh generators allow us to create a pure quadrilateral mesh or to distinguish separately both the areas covered by the most popular triangular elements and the infinite part covered by quadrilateral 8-node elements. Moreover, the authors are not familiar with the generator which would create 5-node infinite quadrilateral mapped elements. Therefore, a self-made generator was used to build a mesh containing only quadrilateral elements.
Fortunately, in optical mammography light sensors and sources are located in a special hemispherical or cone shape constant form, so the effort related to the manual creation of the infinite element mesh, which surrounds the area of interest, has to be made only once.

Damped walls
Measured potential values were compared to these calculated with BEM (used to solve forward problem) and additionally with calculation results taken from FEM.
Inverse problem solution required 7 iterations to find interpolation coefficients related to investigated conductivity function (16).
Potentials -measured and calculatedfor selected projection angles: 7 (electrodes 7-20) and 10 (electrodes10-17), as defined in figure (13), are presented in figures 14 and 15. Projection 7-20 presents the best achieved results and projection 10-17 the worst. Application of Functionally Graded Materials theory into Boundary Element Method used to find spatial distribution of humidity and salinity into damped walls seems to be better idea than using multiregional BEM due to simpler mesh and boundary conditions. Described method will be used to find thomographic image of dumped wall.
Most of articles discusses objects where FGM changes their properties only in one direction [8,13,19,20].
Damped wall presented above represents the graded material where humidity and salinity varies in at least two directions (corners of the buildings will require 3 dimensional model).
Achieved results are quite promising although are worse than ones received during theoretical calculations based on flat capacitor model with known analytical solution or dumped wall model calculated by BEM and FEM where material properties varied only in one direction.
Possible reason of low accuracy in some projection angles (see Fig. 15) can be caused by the insufficient data associated with too few angle projections. Possessed measurements are more suitable for the model where conductivity distribution varies only in one direction corresponding to the thickness of the wall and not its height.
Next stage will be to calculate theoretical model of the wall where material properties varies in two directions x and y (along wall thickness and height respectively) and compare FEM and BEM results.
We believe that new measurements which will include all possible projections angles allows us to receive results with much better accuracy.
Final model will also include infinite boundary elements to avoid an unknown height of wall humidity and possible related errors. BEM is known to be well suited to handle open boundary models containing infinite elements.