Modeling the shape of boundary using NURBS curves directly in modified boundary integral equations for Laplace’s equation

The main purpose of this study is the modification of the parametric integral equations systems (PIES) method, to include the NURBS curves into its mathematical formulas. We want to create an opportunity for defining the boundary shape using NURBS curves without the classical discretization process. We present an inclusion of the NURBS curves into PIES, and its application in process of modeling and solving the boundary value problems, on Laplace’s equation example. The correctness and accuracy of such modeling is confirmed by the examples. Obtained solutions are compared with these obtained using linear segments and B-spline curves. Additionally, we verify the correctness of modeling and modifying the shape using analytical solutions. The NURBS curves significantly expand the possibilities and increase the accuracy of the boundary shape modeling which is an important aspect of the numerical solving of boundary value problems.


Introduction
For many years researchers have used boundary integral equations (BIE) to solve boundary value problems (Banerjee and Butterfield 1981;Bonnet 1995 Eugeniusz Zieniuk ezieniuk@ii.uwb.edu.pl In this paper, we consider an application of NURBS curves not during numerical solution of BIE (like in isogeometric modeling), but in analytical modification of classical BIE. The modification is strictly connected with inclusion of NURBS curves directly into mathematical formulas of BIE. Such approach allows us to exploit all of advantages of the NURBS curves in the process of solving boundary value problems. As the result of this modification, we obtain the generalized parametric integral equations system (PIES). We also present the method of numerical solution of generalized PIES as well as the results of numerical tests. Finally, we compare results of analytical solutions of some examples with ones obtained by generalized PIES to confirm reliability of PIES solutions. Additionally obtained solutions are compared with these obtained using linear segments and B-spline curves. We present the impact of modeling accuracy on obtained solutions, as well.

NURBS curves in modeling the shape of boundary in boundary value problems
NURBS curve of any degree n is described by control points P i (i = 0, 1, . . . , m), corresponding to particular points weights w i (i = 0, 1, . . . , m) and knots t j ( j = 0, 1, . . . , m + n + 1), and is defined by the following formula (Piegl and Tiller 1997): Additionally, the base function N n i (t) of n degree is defined as normalized B-spline blending function and generally presented by the following recursive formula (Piegl and Tiller 1997): where 0 ≤ i ≤ m, 1 ≤ k ≤ n, 0 0 := 0. Correct shape definition is strictly connected with proper selection of the control points (its amount and location) as well as proper definition of weights and knots. Such way of defining seems to be very sophisticated, however significantly increases the accuracy of modeling and allows to define even complex shapes by one closed NURBS curve only.
For example, in Fig. 1, the impact of weights (in particular control points) on NURBS curve of second degree is presented. The curve is defined by nine control points (where P 0 = P 8 ). As we can see, the value of the weight determines the impact of particular point on the shape of the curve. For lower values of weights, the curve lies far away from the points, while higher ones shift the curve closer to the points. Additionally, NURBS curve allows to model the circle shape almost exactly (the radius of the circle is almost equal in each point of curve, using the weights value equal to 0.707 and the curve of second degree only (Fig. 1). Such way of modeling allows for very easy and quick modification of the curve.
Another, very important part, which determine the shape of NURBS curve, are knots. Knots divide the curve into segments which define the shape. The impact of two values of knot vector on the curve shape (with shapes of corresponding base functions N ) are presented in Fig. 2. Similarly as in previous example point P 0 = P 8 . Presented curve is Weights of points P , P , P , P : 1 3 5 7 P 0 P 1 P 2 P 3 P 4 P 5 P 6 P 7 defined by the following vector of knots: (0, 0, 0, 0.25, 0.25, a, b, 0.75, 0.75, 1, 1, 1), where knots described as a and b are varied to show their impact on the curve.
First case a = b = 0.5 ( Fig. 2a) presents the multiple knots, where the curve intersects particular control points. The number of intersected control points and how many times the value of knot should be repeated is strictly connected with the degree of the base functions. For example, in case of second degree curve, we should use doubled knots. Exceptions are knots at the start and at the end of the curve, which ensure closure of the curve. In this case, the number of the same values of knots are increased by one.
The next case a = 0.47, b = 0.53 (Fig. 2b) presents the division of the curve section (between points P 2 and P 6 ) into three subsections (0. 25, 0.47), (0.47, 0.53) and (0.53, 0.75). In this case, the curve no longer intersects the point P 4 . In the next case (0.36, 0.5) (Fig. 2c), we can notice unsymmetrical shape of this curve section as well as of its base functions. The last case a = 0.36, b = 0.64 (Fig. 2d) is defined similarly to the second one, but with increased central subsection. We should note, that knots are another tool which makes modeling and modification of the curve much easier.
As presented in Fig. 2, the boundary was modeled by continuous closed NURBS curve, defined in parametric space with parameter t (0 ≤ t ≤ 1). The endpoints of each segment are defined by knots. Multiple knots have an impact on continuity of the curve in points of segments join. Direct application of such way of boundary modeling in BIE results in automatic adaptation of BIE to modified boundary. It is impossible using NURBS curves to define boundary elements (previously obtained in process of boundary discretization). Therefore, the aim of this paper is an analytical application of presented effective way of boundary modeling instead of boundary integral in classical BIE.

The strategy of analytical inclusion of NURBS curves into mathematical formulas of BIE
We consider a problem described by 2-dimensional Laplace's equation defined in domain Ω and bounded by boundary S. BIE for Laplace's equation can be represented by the integral identity (Banerjee and Butterfield 1981;Burczyński and Fedeliński 1991): where the value of u(x) depends on the location of x and is given by following expression: where for smooth boundary α = 0.5. If x ∈ S, then formula (3) is the classical boundary integral equation (BIE). In identity (3), integrand U (x, y) is known as the fundamental solution for Laplace's equation (Burczyński and Fedeliński 1991): and is represented by the following formula (Burczyński and Fedeliński 1991;Brebbia et al. 1984): (3) is the classical singular solution represented by traction kernels written as follows: The accuracy of the solutions in domain Ω depends on two factors: (1) approximation accuracy of the solutions on the boundary and (2) approximation accuracy of the boundary shape. Hence, it seems reasonable to modify the traditional BIE to separate (in process of numerical solution) the approximation of boundary from the approximation of unknown functions. In other words, it allows to choose the most convenient methods of boundary modeling considering its complexity without any intrusion into approximation of boundary functions and vice versa. Due to the above given reasons, at first, we made analytical modification of BIE to include the shape of the boundary into mathematical formulas of BIE.

Analytical modification of classical boundary integral equations (BIE)
To obtain modification of generalized integral identity (3) for 2D potential problems modeled by Laplace's equation (Zieniuk 2003a, b), we need to apply the Fourier transform. As the result, we obtained: where U (iξ , y) is the transform of the fundamental solution, whereas P k (iξ , y) is the transform obtained in the domain of the Fourier transform, based on fundamental solution. For Laplace's equation, these transforms are presented in the following form: where After substitution (9) and (10) into (8), we obtained the following expression of Fourier transform in the domain: where boundary functionsp(ξ ) andũñ k (ξ ), k = 1, 2 are represented by means of the following boundary integrals: uñ k (ξ ) = S u(y)n k (y)e −i{ξ 1 y 1 +ξ 2 y 2 } dS(y).
It should noted, that in Eq. (12), the boundary is separated from the domain. The boundary is defined in the mathematical formulas of boundary integrals (13) and (14). In our further considerations integral (14) is used to describe transformũñ k (ξ ) on the boundary S. The unknown identity function u(y) in (14) can be defined by means of the following Fourier transform: where integrand u(iω) is written as: The formula (16) is particular case of the transform (12) in a smooth boundary.

The definition of the smooth contour integral by segments
After insertion (15) into (14), we obtained the convolution integral equation in the domain of the Fourier transform:ũñ where the kernel is: (16), and ξ = ω. The contour integral in (18) takes into consideration the shape of boundary S. For further transformation of the above analytical formulas, we need to specify the boundary S. We assume that the boundary is defined by a finite number of the curved segments S l (l = 1, 2, . . . , n) with the necessary continuity conditions in points of segments join. During the process of analytical modification of BIE, we defined these segments using different kinds of curves widely used in computer graphics (Goshtasby 2005;Mortenson 1985Mortenson , 1999. The NURBS curves can be also applied, like in isogeometric analysis (Piegl and Tiller 1997), but at the level of analytical modification of the identity.
Therefore, considering the boundary geometry (defined by segments) in kernel (18), we obtain: where The boundary transformsũñ k (ξ ) on the left side of (17) may be represented in the following form:ũñ whereũñ (l) We define the boundary transformsp(ω) andũñ k (ω) on individual segments on the right side of (16):p After substituting (19), (21) and (23) in (17), we obtain the following system of the convolution integral equations: where l = 1, 2, 3, . . . , n, and: Next, we assume that segments S l (t), which define contours of integrals, are mathematically defined by parametric functions l (t)] T with parameter t. Hence, these functions must be considered in boundary transforms: (20), (22) and (23), in which the boundary is still defined in a general way by means of segment contour integrals S l (t).
The kernel (20), for the boundary defined by parametric segments, is described by the following formula ( and integral (22) by the following formula: where t l−1 , t l are unique values contained in the knot vector and they mean the beginning and the end of the segment with index l, respectively. The integrand functions: were obtained after the substitution of curve segments S (i)

Generalized parametric integral equations system (PIES)
After application of the inverse Fourier transform and proper transformations into (24), we obtained the following formula, which requires further integration: After integrating both sides of Eq. (29) with respect to ξ , we obtained: Next, to calculate integral with respect to t, we inserted S (1) Finally, after substitution and integration, the formula (30) is presented as follows: where To remove the Dirac distribution δ(λ 2 ) from (31), we calculated integrals on both sides of the equation with respect to x 2 . After substitution Φ(x 1 ) = t, Eq. (31) finally is defined as: where l = 1, 2, 3, . . . , n. Further integration of right side of (32) with respect to ω requires definition of integrand function u ( j) (iω) in explicit form. This function is presented by Eq. (25), wherep ( j) (ω) and uñ ( j) k (ω) are defined on segments j similarly like in case of l segments in (27): where s j−1 , s j are unique values contained in the knot vector and they mean the beginning and the end of the segment with index j, respectively, while is the normal vector. The value J j (s) is the Jacobian determined by the following formula: Normal vectors of segments should be directed outside the domain. Declaration of the segments is reduced to defining control points, their weights and knots only. The way of defining these points has an impact on the direction of the normal vector. After substitution (33) and (34) into (25) and next (25) into (32), we obtained the formula named parametric integral equation system (PIES) (Zieniuk 2003a, b) in the following form: where J j is Jacobian (35) and the kernels in formula (36) are represented by functions U * l j and P * l j . The first integrand U * l j is represented by the following expression: where After calculation of integral (37), we obtained an expression in an explicit form, which include the shape of boundary into its formulas (by function η(t, s) of two parameters t and s): As a result of derivations, there were attempts to present (39) in the same form as classical fundamental solution. Nevertheless, in Eq. (39), η 1 and η 2 are not simply variables like in fundamental solution, but they are functions of parameters t and s, like it is presented in (38). Therefore, Eq. (39), contrary to fundamental solution, include the shape of boundary into its formulas and was called as a fundamental boundary solution.
In PIES, indexes l, j = 1, 2, 3, . . . , n allow to define relationship between boundary segments of considered boundary value problems. In this paper, we define the shape of boundary in PIES by NURBS segments. Relationship between segments, that model closed boundary, was presented by integral equations system. However, PIES was not defined on boundary like classical BIE, but in the domain of boundary definition (means 0 ≤ t, s ≤ 1).
An important advantage of such modeling is that for any shape of the boundary modification (as is shown in Sect. 7) it will be defined always in the same domain. This is the domain of NURBS curve definition, which model closed boundary. In other words, PIES is automatically adapted to modified shape of boundary without any additional computations. The modification is connected with the change of control points, their weights and knots only, like in Sect. 2. These changes result in automatic adaptation of PIES to new shape of boundary, because NURBS curves are directly included into mathematical formulas of PIES.
Solution of PIES is reduced to finding boundary functions defined on NURBS segments (included into PIES kernels). Development of global algorithm of solving PIES is much easier, because PIES is not defined on boundary contrary to BIE. PIES is defined in parametric reference system corresponding to domain of NURBS curves, which are used to define the boundary.

Modification of integral identity for solutions in domain
After PIES solution, we obtain only solutions on boundary. To obtain solutions in the domain, proper adaptation of (12) is required, where boundary integralsp ( j) (ξ ) andũñ ( j) k (ξ ) are represented similarly to (23). Finally Eq. (12) is presented as follows: and s j−1 , s j are unique values contained in the knot vector. They mean the beginning and the end of the segment with index j, respectively. After the inverse Fourier transform of (42), we obtained the following expression: The first integrand U * j (x, s) takes the following form: where (45) is represented by the following expression:

Numerical solution of PIES
Mathematical formulas of PIES is defined in the parametric domain of parametric segments. They theoretically described the shape of the boundary. The separation of approximation of the shape of boundary from the approximation of the boundary functions makes possibility to consider these problems independently. To solve PIES, we can use one of methods dedicated to solving the classical singular integral equations. We should remember that first we must define the boundary in PIES. Solution of PIES (36) is reduced to approximation of the boundary functions p ( j) (s) and u ( j) (s) defined on following segments with index j. The given boundary conditions as well as searched boundary functions are approximated on each segment by the following series:

The approximation of PIES
After substituting (48) into (36), we obtained approximated expression of PIES for smooth boundary (where α = 0.5) in the following form: where integrands U * l j (t, s), P * l j (t, s) are represented by formulas (39) and (41). To solve numerically formula (49), pseudospectral method is used. After applying (49) for all collocation points t (c) , located in the parametric segments t l−1 < t (c) < t l of S l (l = 1, 2, 3, . . . , n), we obtained the system of algebraic equations in the following form: The size of the system depends on the number of collocation points. After solving (50), we obtained values of coefficients of terms (48), represented by vectors [h l j ] or [g l j ] (according to posed boundary conditions).
The elements of matrices H and G are calculated on the basis of the following expressions: The matrix elements in (50) requires the calculation of strongly and weak singular integrals. After taking into account the boundary conditions in (50) and solving the algebraic system of equations, we obtained the unknown coefficients u ( j) r or p ( j) r on individual segment of the boundary.

Approximation of integral identity for solutions in domain
After solving equations (50), we obtained solutions on the boundary. Solutions in the domain Ω at point x = (x 1 , x 2 ) are computed by integral identity (45). After substituting approximating series (48) into (45), we obtained the identity in the following form: The integrands U * j (x, s), P * j (x, s) are represented by (46) and (47). To determine the solution in the domain, we need to take into account in the formula (53) only coefficients u

Testing examples
The developed strategy of generalized PIES should be tested in respect of modeling and solving boundary value problem accuracy. First, we decided to compare the solutions obtained using NURBS with solutions obtained formerly in PIES (for linear segments and B-spline curves). We also verified the efficiency of boundary shape modification using NURBS curves. For this purpose, we were changing the location of control points, their weights and knots of the NURBS curve. Obtained solutions were compared with analytical ones. The mentioned tests, on examples of different boundary value problems, were presented below.

Comparison of the PIES solutions between application of linear segments and NURBS curves
First, to verify the correctness of boundary shape modeling in PIES by NURBS curves, we considered example of rectangular domain. modeling of such shape by linear segments as well as by NURBS curves should be made with the same accuracy. In Fig. 3, we present the boundary modeling of such domain by NURBS curves of first, second and third degree. Additionally, we present corresponding normalized B-spline blending functions N . The weights, in each case (Fig. 3a-c), are equal to one. The vector of knots, for the curve of first degree, is defined as: (0, 0, 0.25, 0.5, 0.75, 1, 1), for the curve of second degree: (0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1, 1) and for the curve of third degree (0, 0, 0, 0, 0.25, 0.25, 0.25, 0.5, 0.5, 0.5, 0.75, 0.75, 0.75 , 1, 1, 1, 1). Then, we consider the boundary value problem with the boundary conditions presented in Fig. 4. We present not only control points, but also points in domain in which the solutions were obtained.
In Table 1, we presented the solutions of above example (Fig. 4). Solutions are obtained for different degree of NURBS curves, as shown in Fig. 3 and for linear segments (points are defined like in Fig. 4) (Zieniuk 2002). The accuracy of the linear segments as well as the NURBS segments modeling in this case should be the same. Therefore, obtained solutions (presented in Table 1) are almost the same.

Comparison of PIES solutions between application of B-spline and NURBS curves
Next, we check the impact of curvilinear shape modeling accuracy on accuracy of PIES solutions. For this purpose, we considered curvilinear example presented in Hromadka (1987) based on Laplace's equation. In this example, on elliptical boundary x 2 1 /a 2 + x 2 2 /b 2 = 1, the following Dirichlet boundary conditions were defined: Analytical solutions, for so-defined problem, are presented in Hromadka (1987) in the following form: The way of the boundary shape modeling, using classical B-spline curves and NURBS curves, was presented in Fig. 5. Previously in PIES, to define curvilinear shape of boundary Bspline curves were used (Zieniuk 2007). These curves, as it is presented in Fig. 5a, are defined by the de Boor's approximation control points D i (i = 0, 1, . . . 7). However, in practice, such shape of boundary modeling by approximating points is quite difficult. Therefore, only points located on boundary are defined P i (i = 0, 1, . . . 7) and then de Boor's control points were calculated D i (i = 0, 1, . . . 7).
In case of NURBS curves (Fig. 5b), to model the boundary shape, we should define the points P i (i = 0, 1, . . . 7), where only P i (i = 0, 2, 4, 6) are located on boundary. In addition, weights w i (i = 0, 1, . . . 7) for each point and knots vector (0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4) are defined. In this example, normalized B-spline blending functions are the same as those presented in Fig. 3b. The only difference is the location of control points. Additionally, in   case of points P 1 , P 3 , P 5 and P 7 we define weights w 1 = w 3 = w 5 = w 7 = 0.707 to improve the accuracy of the elliptical boundary modeling. In B-spline application (Zieniuk 2007) the third degree curves were used, while in NURBS the second degree curve only. PIES solutions for both ways of modeling with BEM and analytical solutions were presented in Table 2. Additionally, in Table 2, we present percentage relative error of numerical solutions in domain with respect to analytical solutions.
Average relative error of PIES solutions, in comparison to analytical ones, is much smaller in case of NURBS curves application than in B-spline curves. BEM often used linear shape functions and it is assumed, that for the engineering purposes obtained solutions are satisfactory. However, based on presented tests, we should noted that the accuracy of boundary shape modeling has an significant impact on the accuracy of obtained solutions. Therefore, it is reasonable to use contemporary methods of the boundary shape modeling.

Modification of the boundary using control points and weights
Next, the authors were focused on the opportunities connected with the inclusion of NURBS curves into mathematical formulas of PIES, i.e., easy and convenient way of modeling as well as modifying the shape of boundary in PIES. The mathematical formulas of PIES automatically adapts to the changes of input values such as control points, their weights t t , 9-11 0-2 t 3-4 t 5-6 t 7-8 x 1 x 2 10 (a) (b) Fig. 6 Examples of modeling and modifying the shape of boundary using NURBS curves and knots. First, to illustrate easiness of the boundary shape modification, we considered examples presented in Fig. 6. Presented examples are obtained as a result of elliptic shape modifications (Fig. 5) using control points and their weights. Therefore, the shape is still defined as closed curve of second degree. We decided to change only control points and their weights. The knot vector remains the same. Therefore, functions N are presented similarly to ones in Fig. 3b. In Fig. 6 we marked knots and only these weights, which values are different than one. In Fig. 6a the weights of points P 1 and P 7 are equal to 1. In the next example, presented in Fig. 6b, to shift the curve closer to points P 1 and P 3 , their weights were increased to 3. The other modifications, necessary to obtain these shapes presented in Fig. 6 from ellipse (Fig. 5b), assumed changes of coordinates of control points only.
To verify the accuracy of more sophisticated shapes modeling, we decided to compare PIES solutions of boundary value problem with analytical ones. Boundary shape definitions of the examples were presented in Fig. 6. The boundary value problems were modeled by Laplace's equation. We used the function: which satisfies Laplace's equation. Based on this function we can obtain the Dirichlet and Neumann boundary conditions for any considered shape of boundary. In examples presented in Fig. 6, we considered the Dirichlet boundary conditions u(x). Solutions in the domain were obtained in the same cross-section, where x 2 = 3.5 and x 1 was changed from 1.5 to 6.5. The points of cross-section are presented in Fig. 6. The analytical solutions (u a ) and PIES ones (Fig. 6a, b) are presented in Table 3.

Modification of the boundary using control points, weights and knots
Solutions of problems presented in Fig. 6 were satisfactory. Therefore, we decided to tests examples, where all of NURBS curve input data (control points, their weights and knots vector) were changed. Such modification (based on ellipse, as well) and corresponding N functions are presented in Fig. 7.
Both shapes are defined by the same weights values except for points P 1 and P 5 . To shift the curve closer to points P 1 and P 5 we increased their weights to value 3 (Fig. 7a) and to value 5 (Fig. 7b). Next difference between shapes presented in Fig. 7a, b is the location of P 1 and P 5 control points and changes of some values in knots vector. Double knots t 3−4 and t 5−6 in Fig. 7a have been changed in Fig. 7b to single knots t 3 , t 4 , t 5 and t 6 . The impact of Table 3 Comparison between analytical (u a ) and PIES solutions of problems presented in Fig. 6 x u a Fig. 6a Error (%) Fig. 6b Error (%)   Fig. 7. We should noted, that N functions for other values of knots (t 7−11 ) remained unchanged. In next example (Fig. 8a) we changed the weight of the point P 1 into value 3. Other weights did not modify the shape, therefore all of them were equal to 1. To smooth the curve in points P 2 , P 4 and P 6 , three double knots were split into six single knots. Therefore, the curve does not intersect points P 2 , P 4 and P 6 , contrary to the case of double knots. In the example presented in Fig. 8b the weight of P 1 was also increased (w 1 = 3), while the weight of P 6 remained unchanged. Other weights did not modify the shape, therefore all of them were equal to 1, as well. This example presents splitting of two double knots in points P 2 and P 6 . Therefore, the curve does not intersect these two points. This results in a smooth connection between segments.
Above mentioned examples clearly present, that NURBS curve is an easy tool for defining linear as well as curvilinear segments of the boundary shape using only one closed curve. In addition, changes of values in knots vector are shown as the chart of N functions (respectively to shapes in Figs. 7 and 8). To verify the modeling as well as implementation of NURBS curves in PIES, we solved boundary value problems presented in Figs. 7 and 8. The Dirichlet boundary conditions were defined by (56). Obtained solutions in domain, in cross-section shown in Figs. 7 and 8 as dotted line, are presented in Table 4.
In case of all examples presented in Figs. 6, 7 and 8 the average relative error of PIES solutions in domain, with respect to analytical solutions, is very small. We should noted, that numerical solutions of PIES almost match the analytical solutions, even in the case of more sophisticated shapes of the boundary. Therefore, application of NURBS curves into mathematical formulas of PIES, results in high accuracy of solutions as well as more accurate shape of boundary modeling. In addition, it seems to be very easy and functional tool for modeling and modification of the boundary shape.

Conclusions
In this paper, as result of NURBS curves insertion (widely applied for the shape of boundary modeling) directly into mathematical formulas of classical BIE, the generalized PIES were obtained. It is a different way comparing to known from literature application of NURBS at the stage of numerical solution of BIE (called as isogeometric modeling in BEM). Application of NURBS curves in mathematical formulas of generalized PIES have a great advantage over direct application of curves to define boundary elements. The changes of the shape of Table 4 Comparison between analytical (u a ) and PIES solutions of problems presented in Figs. 7 and 8 x u a Fig. 7a Error (%) Fig. 7b Error (%) Fig. 8a Error (%) Fig. 8b Error (%) boundary in PIES, using control points, their weights and knots, result in automatic adaptation of PIES to modified shape of boundary. Proposed application of NURBS curves in PIES allows us to define the shape of boundary more accurately. We have noted a significant increase in accuracy of obtained solutions comparing to ones previously obtained by application of B-spline curves in PIES. Additionally, for modeling ellipse we can use the second degree NURBS curve only, while in case of Bspline curves cubic segments should be used. Next advantage is the opportunity of modeling of continuous shape of boundary using only one closed curve. Therefore, the application of NURBS in PIES increased the accuracy of modeling of the shape of boundary (it also increase the accuracy of the solutions), as well as creates opportunity of easy and uniform way of defining the shape of boundary.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.