Exact 3D Greens Function Integrations on Triangles

This demo computes exact integrations of 3D Green's functions and compares them with numerical solutions.

Table of Contents

• Introduction
• How to Use
• Demo
• Results
• Discussion
• Acknowledgement
• Citation
• License

Introduction

We compute analytically the Green's and its gradient and product with linear functions on a triangle T in R3.

Let x be an evaluation point, r(x,y) = |y-x| and Let h(x) be a linear function on the triangle with value 1 at the first triangle vertex and 0 on the other vertices.

We compute the following integrals

$\int_T \ G(x,y) dA_y = \int_T \ \frac{1}{|y-x|} dA_y$

$\int_T \ h(x) \ G(x,y) dA_y = \int_T \ h(x) \ \frac{1}{|y-x|} dA_y$

$\int_T \ \partial_{n_y} G(x,y) dA_y = \int_T \ \langle n_y , \nabla \frac{1}{|y-x|} \rangle dA_y$

$\int_T \ h(x) \ \partial_{n_y} G(x,y) dA_y = \int_T \ h(x) \ \langle n_y , \nabla \frac{1}{|y-x|} \rangle dA_y$

$\int_T \ \nabla G(x,y) dA_y = \int_T \ \nabla \frac{1}{|y-x|} dA_y$

$\int_T \ h(x) \ \nabla G(x,y) dA_y = \int_T \ h(x) \ \nabla \frac{1}{|y-x|} dA_y$

The analytic computation is performed based on [Graglia]'s results.

The following reference was the base for the formulae in this code: [Graglia] R. D. Graglia, "Numerical Integration of the Linear Shape Functions Times the 3-D Green's Function or Its Gradient on a Plane Triangle," IEEE Transactions on Antennas and Propagation, vol. 41, no. 10, pp. 1448–1455, Oct. 1993.

How to Use

Copy greens_triangle_integration.py into your project and import it using

import greens_triangle_integration

Function:

greens_triangle_integration.integrate(
    evaluation_points, triangle, values_of_interest)

• Input:
  • evaluation_points
    • A numpy array of the $x \in \mathbb{R}^3$ values.
  • triangle
    • a numpy array of 3 points in $\mathbb{R}^3$ that represent a non-degenrate triangle and defines the integration domain $y\in T$.
  • values_of_interest

    • a list of strings specifying what to return in the dictionary of results.

      Value of Interest	Integral
      'G'	$\int_T \ G(x,y) dA_y$
      'h_G'	$\int_T \ h \ G(x,y) dA_y$
      'n_grad_G'	$\int_T \ \partial_{n_y} G(x,y) dA_y$
      'h_n_grad_G'	$\int_T \ h(x) \ \partial_{n_y} G(x,y) dA_y$
      'grad_G'	$\int_T \ \nabla G(x,y) dA_y$
      'h_grad_G'	$\int_T \ h(x) \ \nabla G(x,y) dA_y$

    • any call involving $h$, returns a numpy array of dimensions (nr_evaluation_points,3) with the $i$-th column being the integral with the linear function $hi$, where $hi$ is the linear function that reaches the value 1 on the $i$-th triangle vertex and 0 on the other vertices.

• Output:
  • integration_results
    • The analytic integral values on the evaluation points, in the form of numpy arrays, stored in a dictionary that can be accessed by the strings mentioned in the input.

Example use:

import numpy as np
import greens_triangle_integration

# prep
evaluation_points = np.random.rand(100, 3)
triangle = np.random.rand(3, 3)
values_of_interest = ['h_G','h_n_grad_G']
# do
integration_results = greens_triangle_integration.integrate(
    evaluation_points, triangle, values_of_interest)
# read 
Int_G = integration_results['h_G']
Int_h_n_grad_G = integration_results['h_n_grad_G']

Demo

Simply run the greens_triangle_integration_demo.py, and pip install missing packages if prompted.

An abitrary triangle is generated by the vertices triangle[0], triangle[1], triangle[2]. A wide selection of test evaluation points are generated. The a numerical integration is performed via scipy.integrate.dblquad.

The test evaluation points are lines going through:

• a point further outside of the triangle plane
• a point in the triangle plane but not in the triangle or along an extended edge
• a point on the extended edge of a triangle
• a point on the extended edge of a triangle but off the other direction
• a point in the triangle but not on an edge
• a point on an edge but not a vertex
• a triangle vertex (the first one, where h=1)

The formulas have been adjusted to not require a transformation into the local triangle frame. Additionally, an orientation dependant singularity issue on the extended edge of a triangle has been fixed.

Important note: When an evaluation point is at the boundary of the triangle, the tangential part of the gradient integrals become singular.

Results

Running the code will show you the following results:

Figure 1: Integration results for a line through a point outside of the triangle plane.

Figure 2: Integration results for a line through a point in the triangle plane but not in the triangle or along an extended edge.

Figure 3: Integration results for a line through a point on the extended edge of a triangle.

Figure 4: Integration results for a line through a point on the extended edge of a triangle but off the other direction. [Graglia]'s formulae would fail in this case.

Figure 5: Integration results for a line through a point in the triangle but not on an edge.

Figure 6: Integration results for a line through a point on an edge but not a vertex. There is a singularity.

Figure 7: Integration results for a line through a triangle vertex (the first one, where h=1). There is a singularity.

Discussion

The analtic integration results are robust to the many possible cases and return the same result as the numerical integration based on scipy's dblquad. Away from singularities the analytic integration is around 10e5 times faster, while closer to the singularities this reaches 10e6 speed-up using the default settings of dblquad.

However, the singularity at the tangential component when integrating $\nabla G(x,y)$ when $x \in \partial T$ seems to be unresolvable. Luckily, the normal gradient $\partial_n G(x,y)$ does not have this problem. The main problem is that the values of f2 in the code become ill defined. At the tangential singularities they are set to zero, which is useful when having two triangle in the same plane next to each other, where these singularities would cancel out.

This code was created for dealing with boundary integral equations.

It is straight forward to apply this code to flat polygons for the integrations not involving h(x). For h(x) integrations the flat polygon must be split into triangles.

Acknowledgement

I thank Fabio Freschi for having uploaded a limited (without h(x)) matlab implementation that helped me get started. int_green3d (https://www.mathworks.com/matlabcentral/fileexchange/47782-int_green3d), MATLAB Central File Exchange. Retrieved February 5, 2025.

Citation

If you use this code in your research, please cite it as follows:

@misc{padilla2025exact3Dgreens,
    author = {Marcel Padilla},
    title = {Exact 3D Greens Function Integrations on Triangles},
    year = {2025},
    publisher = {GitHub},
    journal = {GitHub repository},
    howpublished = {\url{https://github.com/marcelpadilla/exact-3D-greens-function-integrations-on-triangles}},
}

License

This project is licensed under the MIT License. See the LICENSE file for details. Enjoy using it for whatever projects you have. Just please cite it.