{'fraction_non_alphanumeric': 0.0491326895565093, 'fraction_numerical': 0.02405221745350501, 'mean_word_length': 4.121108058608058, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 9, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-gravity in curved mesh elements

The local character of self-gravity along with the number of spatial dimensions are critical issues when computing the potential and forces inside massive systems like stars and disks. This appears from the discretisation scale where each cell of the numerical grid is a self-interacting body in itself. There is apparently no closed-form expression yet giving the potential of a three-dimensional homogeneous cylindrical or spherical cell, in contrast with the Cartesian case. By using Green's theorem, we show that the potential integral for such polar-type 3D sectors -- initially, a volume integral with singular kernel -- can be converted into a regular line-integral running over the lateral contour, thereby generalising a formula already known under axial symmetry. It therefore is a step towards the obtention of another potential/density pair. The new kernel is a finite function of the cell's shape (with the simplest form in cylindrical geometry), and mixes incomplete elliptic integrals, inverse trigonometric and hyperbolic functions. The contour integral is easy to compute; it is valid in the whole physical space, exterior and interior to the sector itself and works in fact for a wide variety of shapes of astrophysical interest (e.g. sectors of tori or flared discs). This result is suited to easily providing reference solutions, and to reconstructing potential and forces in inhomogeneous systems by superposition. The contour integrals for the 3 components of the acceleration vector are explicitely given.

Consequently, in continuous media, the integral approach cannot provide reliable potential values inside matter from simple quadrature schemes unless kernel singularities are properly treated, either by direct integration or through specific techniques. While the Poisson equation appears as the privileged alternative, it requires the knowledge of accurate boundary conditions that only the integral approach can provide, except in the very specific case of periodic boundary conditions. The gradual increase in the precision of models and simulations is a considerable source of motivation to derive new formulae, and explore various techniques and algorithms to solve this problem more and more efficiently (Grandclément et al, 2001;Matsumoto and Hanawa, 2003;Huré, 2005;Jusélius and Sundholm, 2007;Li et al, 2008;Guillet and Teyssier, 2011;Ma et al, 2012).
This article deals with the computation of the gravitational potential inside continuous systems, homogeneous or not, that can be discretised on a polar-type 3D-grid made of homogeneous cells/sectors. This includes cylindrical and spherical grids which are commonly used in numerical models and simulations of bodies deformed by rotation and which exhibit a curved shape, like stars and disks. Applications however exceed the astrophysical context. More precisely, we show that the potential of an elementary sector of such grids is always reducible to an onedimensional, contour integral. Dafa-Alla et al (2010) have studied this issue from a finite difference approach. Our work extends a result already known under axial symmetry for closed toroids (see e.g. Ansorg et al, 2003). Unfortunately, we failed to derive an algebraic, close-form expression for this contour integral. However, the integral formula is relatively simple and general, and can serve to generating reference solutions or accurate potential values at some points of space. From a purely practical point of view, the present formula is superior (accuracy and computing time) to the multipole approach which necessarily requires one triple integral over the material volume per term in the infinite series.
The paper is organised as follows. In Section 2, we recall Newton's integral formula for curved, polar-type domains and set the notations and hypothesis. In Section 3, we discuss the conversion of the potential integral into a one-dimensional, contour integral through Green's theorem. We briefly comment on the properties of the new kernel. We then apply the result to the cylindrical cell; this is the aim of Section 4. Explicit formulae as well as several examples are given. An appendix contains a few formulae and a demonstration. A basic Fortran 90 program is available as a supplementary resource linked to the online version of the paper (see Appendix C). Specific computations are also possible upon request. Expressions for the 3 components of the acceleration are given in the Appendix D.

Theoretical background
A volume V of space containing a mass density ρ generates at a point P(r) a scalar potential (Kellogg, 1929): where r ′ refers to points P ′ belonging the source, d 3 τ is the elementary volume, G for Gravitation problems. Successive integrations of the Green function 1/|r − r ′ | over the three spatial directions can sometimes lead to an algebraic formula (e.g. MacMillan, 1930;Durand, 1953;Binney and Tremaine, 1987), but this is rare. Then the Green function, which hyperbolically diverges inside sources, is usually expanded into spherical harmonics. In cylindrical coordinates where P(R, θ, Z), P ′ (a, θ ′ , z) and d 3 τ = adadθ ′ dz, this expansion takes the classical form (e.g. Cohl and Tohline, 1999): where Jm is the Bessel function and m is integer. In the absence of any symmetry, the potential is therefore given, at any point of space, by an infinite series where each term is a quadruple integral. From a computational point of view, this is obviously prohibitive. Besides, the truncation of the series also generates an error that adds to that due to the quadrature scheme. Many terms in the series -tens or hundreds-are necessary (Hachisu, 1986;Stone and Norman, 1992;Mach and Malec, 2012). As often quoted in the literature, the oscillatory character of Bessel functions together with their range of definition are not well suited for numerical calculus. For all these reasons, we easily understand that such an expansion is not especially recommended, even if there is often no way around it. This remark is not specific to Eq.(2) and holds for spherical harmonics as well.
In the case of homogeneous systems as considered here, the number of integrals can advantageously be lowered. As a matter of fact, the integration over the polar angle θ ′ in Eq.(1) is not trivial but feasible (Durand, 1953;Zhu, 2005), which renders for instance Eq.(2) obsolete in the case where the system is θ ′ -invariant (i.e., ρ does not depend on θ ′ inside V whatever a and z). Actually, we can rewrite Eq.(1) as: is a quantity in the range [0, 1]. Assuming that the volume V is a curve sector bounded by the two polar angles θ ′ 1 and θ ′ 2 , then the above expression can be rewritten in terms of incomplete elliptic integrals of the first kind by setting: So, we get: where is the incomplete elliptic integral of the first kind (Gradshteyn and Ryzhik, 2007;Olver et al, 2010), k is the modulus defined by Eq.(4), S denotes the lateral face of the material volume, 2β 1 = π − (θ ′ 1 − θ) and 2β 2 = π − (θ ′ 2 − θ) according to Eq.(5).  We have now a bi-dimensional integral. Equation (6) applies to a wide variety of shapes, including cylindrical and spherical cells, provided the mass density is uniform along any arc. An example of four possible shapes of astrophysical interest is shown in Fig. 1. Note that the integrand in Eq.(6) is basically the potential of a circular arc (Zhu, 2005). The computation of F (β, k) must be performed with caution as soon as the amplitude β stands outside the range [0, π 2 ] (see Appendix A). In the axially symmetrical case, we have for instance θ ′ 1 −θ = 0 and θ ′ 2 −θ = 2π, and so, as expected (Durand, 1953), F (β 1 , k) − F (β 2 , k) = 2K(k) where K(k) ≡ F ( π 2 , k) is the complete elliptic integral of the first kind.
It is in principle possible to compute ψ(r) from Eq.(6) for a given geometry by performing a double numerical integration over a and z. Not only this remains computationally costly, but there is also a logarithmic singularity to manage (due to F ) as soon as the amplitude β = π/2 and k = 1. If the body is a continuous pileup of polar cells along the z-direction, then we can use the closed-form for the potential of the polar cell (Huré, 2012). In this case, there is a single numerical integration to perform (see Sect. D).

Green's theorem
It is clear from Eq.(6) that dadz is an area element. There is therefore an alternative means since the surface S is, following our hypothesis, fully planar. Actually, the double integral can be converted into a line integral by using Green's theorem provided the kernel is the normal component of the curl of a vector W, namely: where d 2 A = d 2 A n is an area element oriented outward, and dℓ is an infinitesimal displacement counter-clockwise along the boundary ∂S of S in the plane ⊥ n. In the present situation, the challenge is to find the appropriate vector W with cylindrical components (M, 0, N ) such that: Clearly, M and N depend on a and z. In this case, we have which is therefore a one-dimensional integral (here, ρ 0 is a mass density, introduced for homogeneity reasons). The existence of the two functions M and N is not guaranteed at all, not only due to the presence of the special function F , but also because ρ varies in general inside the cell (except along any arc). In numerical simulations, the elementary cell is precisely the smallest resolved pattern which marks the limit of the numerical discretisation; in this context, the cell is, most of the time, considered as a homogeneous system, i.e. ρ = ρ 0 . This is the assumption we will make in the following: ρ is therefore constant in the whole integration domain V.
As the potential depends on the difference of two incomplete elliptic integrals with particular amplitudes β 1 and β 2 , we see that we can attack the problem by solving the generic problem: where β 0 is some amplitude and the two functions M 0 and N 0 are to be determined. In these conditions, the two functions M and N defined in Eq.(10) are given by: where the subscripts 1 and 2 refer to the two amplitudes β 1 and β 2 , respectively. We could start the resolution of this problem by considering a general expression for M 0 and N 0 on the basis of a linear combination of the incomplete elliptic integrals, namely F (β 0 , k), the incomplete elliptic integral of the second kind E(β 0 , k), and possibly the incomplete elliptic integral of the third kind Π(β 0 , m 2 , k) with characteristic m. However, we are helped in this task by remembering that the problem has already been solved under axial symmetry (Ansorg et al, 2003). Actually, for a maximum opening angle where ζ = Z − z and E(k) = E( π 2 , k). We see the great advantage of this formulation: both M and N are finite everywhere 1 , even onto the boundary ∂S where k reaches unity, making the determination of the potential through Eq.(10) particularly simple and efficient.
We then use this nice and powerful result -reproduced in details in Appendix B -as the starting point of our investigation. By analogy, we build the incomplete form of the axially symmetrical solution, namely: where f and g are unknown functions to be determined (but which must be zero under axial symmetry). The factor 1 2 -not present in Eq.(11) -comes from the fact that M 1 = −M 2 under axial symmetry (see the end of Section 1). The same kind of strategy has been used in (Huré, 2012). This choice is also motivated by the fact that the partial derivatives of the incomplete elliptic integrals are very close to the partial derivatives of the complete ones. Actually, we have (Gradshteyn and Ryzhik, 2007): and and so the differences between the complete version and the incomplete version come uniquely from the last term in Eq.(16) containing trigonometric functions. As we have where (a + R)m = 2 √ aR, and 1 From the asymptotic behaviour of the K-function (Gradshteyn and Ryzhik, 2007), the leading terms are: as k → 1 (corresponding to a → R and z → Z).
we can easily take the partial derivatives of M 0 and N 0 defined above and set constraints upon the functions f and g. After some algebra and rearrangement of terms, we find the general result: In order to satisfy Eq.(11), we conclude that the sum of the last three terms in this relationship must vanish. As we have it turns out that the following two basic pairs (f, g) are appropriate: or any linear combination, for instance of the form αf and (1 − α)g where α ∈ [0, 1] (any function of a can be added to f , and any function of ζ can be added to g).
The optimal choice may depend on the actual problem. Note that f and g are regular (i.e. non-diverging) functions. Actually, f = 0 as soon as i) ζ = 0 (i.e. in a horizontal plane containing the top/bottom of the cell), or ii) β 0 = {0, π 2 } (i.e. along the line θ = θ ′ 0 ). For a = R and β 0 → π 2 , we also have: The components of the generic problem are finally of the form: We have plotted in Fig. 2 (left panels) M 0 versus Z for R = a and θ − θ ′ 0 = {− π 6 , 0, + π 6 }, and M 0 versus R for θ = θ ′ 0 and Z = { z 2 , z, 2z}. This computation includes the function f with α = 1, and the parameters are a = 1 and z = 1. The function N 0 (including the function g with α = 0) is plotted in the same figure (right panels). We see that M 0 and N 0 are fully regular. ; α = 1 in both cases (i.e. f is included). Right panels: same but for N 0 and α = 0 (i.e. g is included).

The potential of a cylindrical cell
From Eqs. (10) and (25) together with the right expressions for f and g, we can construct the potential due to a given cylindrical cell defined by 6 parameters: β 1 , β 2 , a 1 , a 2 , z − and z + (see Fig. 1a). Fortran 90 routines and a driver program which performs the numerical computation are available as a supplementary resource (see Appendix C). To remove any ambiguity, it is necessary to make explicit the dependency of all the functions with various parameters, i.e. M 0 ≡ M 0 (a, z, β 0 ), N 0 ≡ N 0 (a, z, β 0 ), f ≡ f (a, z, β 0 ) and g ≡ g(a, z, β 0 ). We still omit the position variables R and Z all these functions depend on, as well as the α-parameter (α = 1 to include f and exclude g, and α = 0 in the reverse case). The potential of the homogeneous cylindrical cell is: where and N ≡ N (a, z, β 1 , β 2 ) (28) This formula is a line integral with regular integrands and is exact. As quoted above, this holds even if the top and bottom edges are curved (i.e. if a depends on z), as for instance considered in Ansorg et al (2003). In the present case, ∂S is made of two rectangles, the one located at θ ′ 1 , the other at θ ′ 2 , and each with corners at (a 1 , z − ), (a 2 , z − ), (a 1 , z + ) and (a 1 , z + ), as shown in Fig. 3. To perform the integration along the whole boundary (counter-clockwise), ∂S must be decomposed into four parts: ∂S − and ∂S + for the bottom and top edges (z = z − and z = z + ), and ∂S 1 and ∂S 2 for the left and right edges (a = a 1 and a = a 2 ). We then get the following expression: where we have set: which corresponds to ∂S − and ∂S + , and which corresponds to ∂S 1 and ∂S 2 . There is little chance that these integrations can be performed analytically to produce a closed-form, especially because of the presence of special functions. Figure 4 displays ∆M versus a and ∆N versus z for the following cell parameters: a 1 = 1, a 2 = 3, z − = −1, z + = 1, θ ′ 1 = 0 and θ ′ 2 = π 6 . The computation is performed at the centre of the cell (i.e. R = 2, Z = 0, θ = π 12 ). The area under the curves ∆M (a) and ∆N (z) therefore yields the potential at this particular point.
The numerical integration of ∆M (a) and ∆N (z) is easy to perform and contains no pitfall. Here, we use the second-order trapezoidal rule with n nodes per direction. Figure 5 displays 3 preferential directions crossing-over the cylindrical cell, namely: direction 1, corresponding to the radial direction (R varies while θ and Z are held fixed), -direction 2, corresponding to the azimuthal direction (θ varies while R and Z are held fixed), -direction 3, corresponding to the vertical direction (Z varies while R and θ are held fixed).
The potential ψ(R, θ, Z) has been determined along these 3-directions from Eq.(29) using the basic, trapezoidal rule (the aim is not to optimise the accuracy which can obviously be easily increased by using other schemes). It is shown in Fig.  6. We see that the potential has finite gradients when crossing-over the boundary of the cell, which is typical of matter continuously distributed in three dimensions. The curvature effect is clearly visible from the asymmetrical variation of ψ with R  (direction 1). Left panels in Fig. 7 displays the potential of the cell in the Z = 0, Z = z + and Z = 2z + planes. Right panels are for an opening angle ∆θ ′ = + 4π 3 . Values are given by unit of mass 1 is the opening angle and a 0 = 1 2 (a 2 + a 1 ) is the mean radius. and Z = 2z + (bottom panel). The colour code is the same. The section of the cell is shown in white. The total mass is M ≈ 4.19. Right panels: same but with θ ′ 2 = −θ ′ 1 = + 2π 3 . The total mass is M ≈ 33.5. It happens that there is a closed form -not reported here -for the integration over z in Eq.(6) in the homogeneous case, which corresponds to the potential of a polar cell (Huré, 2012). Thus, the potential of a cylindrical cell can also be obtained numerically by integrating that formula along the z-direction and it is expected to perfectly match the contour integral given by Eq.(26). We can therefore check the correctness of the contour integral relative to direct integration. This is shown in Fig. 8 where we have plotted the relative deviation versus the radius R between the two potential values computed along a quarter of the boundary (i.e., for θ = θ ′ 1 and Z = z − ). We have considered here two different opening angles, ∆θ ′ ∈ { π 6 , 4π 3 }, as for Figure 7. The averaged deviation computed along the boundary versus the number of node n used for the quadratures is shown on the second panel. We see that the two approaches give the same results, the accuracy being limited by the quadrature scheme, which is second-order. In fact, the number of evaluations of special functions is almost the same (although there is no need here for the elliptic integral of the third kind here), and the equation for the boundary z(a) is required in both cases. Finally, acceleration is easily deduced from ψ, and we see from Eq.(26) that the acceleration vector is also defined by a line integral, namely To get the three cylindrical components, we therefore need to calculate the partial derivatives of M and N with respect to R, θ and Z. The final formulae are accessible from the gradients of M 0 and N 0 which are reported in the appendix D. We however note a difficulty here: the acceleration at the edges of the polar cell is logarithmically singular and the formulae given in the appendix have the same drawback. A closed form would really be needed. In the meanwhile, it is always possible to determine the acceleration of the 3D-cell -a finite vector -from the potential through finite differences, which is technically easy.

Conclusion
In this article, we have shown that the gravitational potential of elementary cells which compose certain polar-type grids can be determined from a one dimensional integral. This result generalises an integral expression already known under axial symmetry (Hill and Wheeler, 1953;Ansorg et al, 2003) and enables to treat fully inhomogeneous systems. The new formula is a contour integral which depends on the shape of the grid cells. It means that we can compute it once for all. This is quite advantageous when working in a fixed cylindrical grid. As its kernel is regular, it can be easily computed by standard means with great accuracy. This is especially interesting if one wishes to generate reference solutions for a given cell's geometry. In particular, we have applied the integral formula to the case of a cylindrical cell. The approach can be generalised to other shapes, as soon as they are well bounded in the azimuthal direction. Such objects include the spherical cell, as well as objects with variable vertical thickness as a function of the radius (sectors of tori or flared disks, etc.).
A Values of F (β, k) and E(β, k) for any amplitude β The computation of F (β, k) must be performed with caution as soon as the amplitude β stands outside the range [0, π 2 ]. It is in particular necessary to use the following rules: and even n, and and even n, E(β, k) = nE(k) − E n π 2 − |β|, k , if n π 2 − |β| ∈ [0, π 2 ] and odd n. In order to determine M and N in the following equation we set without loss of generality: where f , g, h and l are four functions to be determined. From Eqs. (18) and (19), and given (Gradshteyn and Ryzhik, 2007): and we get: Forming ∂aN − ∂zM , and gathering terms, we get for the term multiplying kE(k), and for the term multiplying kK(k). The solution by (Ansorg et al, 2003) which corresponds to is obtained for the following settings: which eliminates the term kE(k) and produces the factor 2 a R for the term kK(k).

C A basic Fortran 90 program
Fortran 90 routines and a driver program which computes the M and N functions for a cylindrical cell and the associated potential at one space point (R, θ, Z) are available from the online version of the paper. The quadrature is performed from a Newton-Cotes, secondorder quadrature scheme. External calls to functions IEF(BETA,K) and IEE(BETA,K) refer to the values of the incomplete elliptic integral F (β, k) and E(β, k) respectively which can be obtained from any mathematical library. For single (double) precision computations, change AP into 4 (resp 8), and EPSMACH is the corresponding precision (about 2 × 10 −16 in double precision). These routines are not optimized. Running the code with the default parameters generates the following output: A0,A1,A2,DELTAA 1.0000000 0.75000000 1.2500000 0.50000000 Z0,Z1,Z2,DELTAZ 0.0000000 -0.25000000 0.25000000 0.50000000 THETA0,THETA1,THETA2,DELTATHETA 0.0000000 -0.25000000 0.25000000 0.50000000 MASS 0.12500000 R,THETA,Z,POT 1.0000000 0.0000000 0.0000000 -0.59440136

D Accelerations
Since M and N are build from the two elementary kernels M 0 and N 0 from Eq.(25), we will only expand ∂ R M 0 , ∂ R N 0 , ∂ θ M 0 , etc. After calculus, we get