# Can you calculate the common volume on excel using Boolean Operation

I would like to calculate the common volume of a solid like the one below but at different angles. I can achieve this using ANSYS Design but I was wondering whether it is possible to calculate this on excel or could I program the Boolean Operation into fortran at all? How does the Boolean operation work on ANSYS to calculate the common volume? Thank you for your help!! Essentially, the calculation required is: $$||\Omega|| := \iiint\limits_{\Omega} 1\,\text{d}x\,\text{d}y\,\text{d}z$$ with: \Omega := \left\{ (x,\,y,\,z) \in \mathbb{R}^3 : \begin{aligned} & x^2+y^2+z^2 - \frac{\left(l_1\,x+m_1\,y+n_1\,z\right)^2}{l_1^2+m_1^2+n_1^2} \le R_1^2\,, \\ & x^2+y^2+z^2 - \frac{\left(l_2\,x+m_2\,y+n_2\,z\right)^2}{l_2^2+m_2^2+n_2^2} \le R_2^2\,, \\ & x^2+y^2+z^2 - \frac{\left(l_3\,x+m_3\,y+n_3\,z\right)^2}{l_3^2+m_3^2+n_3^2} \le R_3^2 \end{aligned} \right\} where the integration domain consists of the intersection of three cylinders with axis passing through the origin, of direction $$(l_i,\,m_i,\,n_i) \ne (0,\,0,\,0)$$ and of radius $$R_i > 0$$, with $$i = 1,\,2,\,3$$.

Now, using Wolfram Mathematica 12.0, defining:

{l1, m1, n1, R1} = {1, 0, 0, 1};
{l2, m2, n2, R2} = {0, 1, 0, 1};
{l3, m3, n3, R3} = {0, 0, 1, 1};

cyl1 = (x^2 + y^2 + z^2) - (l1 x + m1 y + n1 z)^2 / (l1^2 + m1^2 + n1^2) - R1^2;
cyl2 = (x^2 + y^2 + z^2) - (l2 x + m2 y + n2 z)^2 / (l2^2 + m2^2 + n2^2) - R2^2;
cyl3 = (x^2 + y^2 + z^2) - (l3 x + m3 y + n3 z)^2 / (l3^2 + m3^2 + n3^2) - R3^2;

Ω = ImplicitRegion[cyl1 <= 0 && cyl2 <= 0 && cyl3 <= 0, {x, y, z}];


and writing:

ContourPlot3D[{cyl1 == 0, cyl2 == 0, cyl3 == 0}, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]
RegionPlot3D[cyl1 <= 0 && cyl2 <= 0 && cyl3 <= 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]


we can get an idea about the layout of the cylinders and the shape of $$\Omega$$: while writing:

Integrate[1, {x, y, z} ∈ Ω]
NIntegrate[1, {x, y, z} ∈ Ω, PrecisionGoal -> 10, WorkingPrecision -> 10]


it's possible to calculate the $$\Omega$$ measure:

8 (2 - √2)

4.686291501

where the exact result can only be obtained in very particular cases such as the one shown in this example. Having understood all this, the rest comes by itself, as it's very simple both to change the angles and radii of each cylinder and to add any other cylinders.

That said, right now I would have no idea how to write an algorithm to be implemented, for example, in VBA in Microsoft Excel, so as to calculate the volume numerically.

This is not an answer which tells you how to compute the volume numerically in any programming language. Instead, I will present some results that may help.

1. For intersection of any number of infinite long cylinders (not necessarily circular), how to express the volume as a sum of line integrals over the edges.
2. In the special case where all cylinders are circular, have same radius and their axes intersect at a single point, how to reduce the sum of line integrals to a sum over the vertices.
3. As a corollary of $$(2)$$, a closed form expression for the case of three cylinders with common radii is derived.

Part I - express volume as a sum of line integrals over the edges.

Let $$p = (x,y,z)$$ stands for any generic point on $$\mathbb{R}^3$$.

Let $$\Omega_1, \Omega_2, \ldots, \Omega_n$$ be $$n$$ infinite long cylinders whose axes are pointing along directions of unit vectors $$a_1,a_2,\ldots, a_n$$. Each $$\Omega_i$$ is invariant under translation along direction of $$a_i$$. We will assume the cross-section of $$\Omega_i$$ cut by any plane perpendicular to $$a_i$$ is a region bounded by a simple closed piecewise $$C^1$$ curve.

Let $$\Omega = \bigcap_{i = 1}^n \Omega_i$$ be their intersection. Its boundary can be composed into a collection of faces $$\partial \Omega = \bigcup_{F\in \mathcal{F}} F$$ and each face lies on one of the cylinders. For any face $$F$$, let $$\Omega_F$$ be the cylinder it lies on (i.e $$F \subset \partial\Omega_F$$). Let $$a_F$$ be corresponding unit vector and $$A_F(p) = (a_F\cdot p) a_F$$. It is easy to verify

$$\nabla \times (p \times A_F(p)) = p - 3 A_F(p)$$

For any $$p$$ on $$F$$, $$a_F$$ is perpendicular to the area element $$dS$$, this means $$A_F(p) \cdot dS = 0$$.

Since $$\nabla\cdot p = 3$$, by divergence theorem and stroke's theorem, we have

\begin{align} V \stackrel{def}{=} \verb/Vol/(\Omega) &= \int_\Omega dV = \frac13\int_\Omega \nabla\cdot p dV = \frac13\int_{\partial \Omega} p \cdot dS = \frac13\sum_{F\in\mathcal{F}}\int_F p \cdot dS\\ &= \frac13\sum_{F\in\mathcal{F}}\int_F (p - 3A_F) \cdot dS = \frac13\sum_{F\in\mathcal{F}}\int_F \nabla\times( p \times A_F) \cdot dS\\ &= \frac13\sum_{F\in\mathcal{F}}\int_{\partial F} (p \times A_F) \cdot dp \end{align} To proceed, decompose each $$\partial F$$ as a collections of edges. Let $$\mathcal{E}$$ be the union of all these edges, we have

$$\bigcup_{e\in \mathcal{E}} e = \bigcup_{F\in \mathcal{F}} \partial F$$

For any edge $$e \in \mathcal{E}$$, it is in contact with two faces in $$\mathcal{F}$$. Pick an arbitrary orientation for $$e$$. Let $$F(e)$$ be the face whose induced orientation on $$\partial F(e)$$ is the same as the one on $$e$$. Let $$G(e)$$ be the other face, the induced orientation on $$\partial G(e)$$ will be opposite to that of $$e$$.

In terms of them, we can express the volume as a sum of line integral over the edges:

$$\bbox[border:1px solid blue;padding: 1em]{ V = \frac13 \sum_{e\in \mathcal{E}} \int_e (p \times (A_{F(e)} - A_{G(e)}))\cdot dp }\tag{*1}$$

For general configuration, this is as far as we can go.

When we restrict ourselves to circular cylinders whose central axes intersect at a single point, the situation is slightly better. The line integrals have the form of incomplete elliptic integral of second kind. However, there are many ways the edges can connect to each other (e.g. they can form loops). A generic simplification seems impossible.

Part II - reduction of line integral to a sum over vertices.

However, there is one special case where improvement is possible. Namely, when all the circular cylinder share the same radius. WOLOG, consider the case where the common radius is $$1$$ and the central axes interest at origin. i.e.

$$\Omega_i = \left\{ p \in \mathbb{R}^3 : |p|^2 - (a_i \cdot p )^2 \right\}\quad\text{ for }\quad 1 \le i \le n$$

Let $$\mathcal{V}$$ be the collection of vertices of $$\Omega$$. For any $$v \in \mathcal{V}$$, Let $$E(v)$$ be the collection of edges connect to $$v$$. For any $$e \in E(v)$$, let $$u$$ be the other end-point of $$e$$. Choose the orientation of $$e$$ so that it is pointing from $$u$$ to $$v$$. Let $$F = F(e)$$, $$G = G(e)$$ to the two faces in contact to $$e$$ and $$a = a_F$$, $$c = a_G$$ be corresponding vectors.

For any $$p \in e \subset \Omega_F \cap \Omega_G$$, we have $$|p|^2 - (a\cdot p)^2 = |p|^2 - (c\cdot p)^2 = 1 \implies (c\cdot p)^2 = (a\cdot p)^2$$ Flipping the sign of $$c$$ is needed, we can assume $$a\cdot p = c \cdot p$$ along $$e$$. i.e. $$e$$ is lying on the plane $$p\cdot( a - c) = 0$$. Choose a coordinate system so that $$a = (0,0,1)$$ and $$c = (\sin\beta,0,\cos\beta)$$.

In this coordinate system, $$\begin{array}{crcl} & |p|^2 - (a\cdot p)^2 = 1 &\iff& x^2 + y^2 = 1\\ \text{ and }\quad & p \cdot a = p \cdot c &\iff& z = \cos\theta z + \sin\beta y \end{array}$$ We can parameterize $$e$$ by a single $$\theta$$. $$p = (x,y,z) = (\cos\theta,\sin\theta,k\sin\theta) \quad\text{ where } \quad k = \frac{\sin\beta}{1 - \cos\beta}$$

Notice \begin{align}(p \times (a - c))\cdot dp &= \left| \begin{matrix} x & y & z\\ 0 &-\sin\beta & 1 - \cos\beta\\ dx & dy & dz \end{matrix} \right| = (1-\cos\beta) \left| \begin{matrix} x & y & ky\\ 0 &-k & 1\\ dx & dy & dz \end{matrix} \right|\\ &= 2(ydx - xdy) = -2d\theta\end{align}

The line integral over $$e$$ can be evaluated as

\begin{align}\int_e (p \times (A_F - A_G))\cdot dp &= \int_e (a\cdot p) (p \times (a - c))\cdot dp\\ &= -2k \int_e \sin\theta d\theta = 2k \int_e d\cos\theta = 2k \int_u^v dx\\ &= \frac{2 c \times a }{1 - c\cdot a}\cdot (v - u) \end{align} Using these, we can decompose the line integrals in $$(*1)$$ into sum of contribution from end points of the edges. Regroup them according to the vertices, we convert the sum over edges to a sum over the vertices.

$$\bbox[border:1px solid blue;padding: 1em]{ V = \frac23 \sum_{v\in \mathcal{V}} v \cdot \sum_{e \in E(v)}\frac{a_{G(e)} \times a_{F(e)}}{1 - a_{F(e)}\cdot a_{G(e)}} }\tag{*2}$$

Part III - closed form expression for volume of three cylinder, same radii.

Back to our problem of three cylinders and consider the special case where all radii are $$1$$.

Let $$a = a_1,b = a_2,c = a_3$$ be the three unit vectors. Define $$\begin{cases} A = b\cdot c\\ B = c\cdot a\\ C = a \cdot b\end{cases} \quad\text{ and }\quad \begin{cases} \tilde{a} = b \times c\\ \tilde{b} = c \times a\\ \tilde{c} = a\times b\end{cases}$$ It is easy to verify $$\begin{cases} |\tilde{a}|^2 = 1 - A^2,\\ |\tilde{b}|^2 = 1 - B^2,\\ |\tilde{c}|^2 = 1 - C^2, \end{cases} \quad\text{ and }\quad \begin{cases} \tilde{a}\cdot\tilde{b} = AB-C\\ \tilde{b}\cdot\tilde{c} = BC-A\\ \tilde{c}\cdot\tilde{a} = CA-B\\ \end{cases}$$ In general, when $$a,b,c$$ are linear independent, $$\Omega$$ will have $$14$$ vertices. $$6$$ of them has degree $$4$$ and the remaining $$8$$ has degree $$3$$.

In order to visualize the configuration, let's imagine $$a,b,c$$ lies close to the +ve $$z$$, $$x$$ and $$y$$ axes. In this case, one of the $$6$$ vertices with degree $$4$$ will lie close to +ve $$x$$-axis. It is $$v_{100} = \frac{\tilde{b}}{|\tilde{b}|}$$. The $$(a_{F(e)},a_{G(e)})$$ for the 4 edges connect to $$v_{100}$$ are $$(a,c), (-c, a), (-a,-c),(c,-a)$$. The contribution of $$v_{100}$$ to $$V$$ in formula $$(*2)$$ equals to

$$\frac43 \frac{c \times a}{|c\times a|} \cdot \left(\frac{c \times a}{1 - c\cdot a} + \frac{a \times -c}{1 + c\cdot a}\right) = \frac83 \frac{1}{\sqrt{1-B^2}}$$

The contribution from other $$5$$ vertices has similar form, with $$B$$ may be replaced by $$A$$ or $$C$$.

An example of a vertex will degree $$3$$ is located at

$$v_{111} = \frac{\tilde{a} + \tilde{b} + \tilde{c}}{\lambda} \quad\text{ where }\quad \lambda^2 = 2(1-A)(1-B)(1-C)$$

The $$(a_{F(e)},a_{G(e)})$$ for the 3 edges connect to it are $$(a,b), (b, c), (c,a)$$. The contribution to $$V$$ from the first edge is

$$\frac23\frac{\tilde{a} + \tilde{b} + \tilde{c}}{\lambda}\cdot \frac{b \times a}{1 - a\cdot b} = -\frac{2}{3\lambda(1-C)}(\tilde{a} + \tilde{b} + \tilde{c})\cdot \tilde{c}\\ = -\frac{2}{3\lambda(1-C)}(AC-B + BC - A + 1 - C^2) = - \frac{2}{3\lambda}(1 + C - A - B)$$ The contribution from other two edges is similar. The total contribution of $$v_{111}$$ to $$v$$ is

$$-\frac{2}{3\lambda}(3 - A-B-C) = -\frac{2}{3}\sum_{(u,v,w) \in \Lambda} \sqrt{\frac{1-u}{2(1-v)(1-w)}}$$ where $$\Lambda$$ is the set of cyclic permutations of $$(A,B,C)$$

The contribution to the other $$7$$ vertices of degree $$3$$ is similar. The end result is

$$\bbox[border: 1px solid blue;padding: 1em]{ V = \frac{4}{3}\sum_{(u,v,w)\in \Lambda'} \left[\frac{1}{\sqrt{1-u^2}} - \sqrt{\frac{1-u}{2(1-v)(1-w)}}\right] }\tag{*3a}$$ where $$\Lambda'$$ is following set of twelve $$3$$-tuples.

\begin{align} &(A,B,C), (-A,-C,B), (A,-B,-C), (-A,C,-B),\\ & (B,C,A), (-B,-A,C), (B,-C,-A), (-B,A,-C),\\ &(C,A,B), (-C,-B,A), (C,-A,-B), (-C,B,-A) \end{align}

In the special case when $$a,b,c$$ are orthogonal to each other, $$A = B = C = 0$$ and the volume reduces to the known result:

$$V = \frac{4}{3} (12)\left[ 1 - \frac{1}{\sqrt{2}}\right] = (16-8\sqrt{2})$$

When $$a, b, c$$ are linear dependent, above arguments cease to work. $$\Omega$$ now has $$2$$ vertices and $$6$$ edges. However, since RHS of $$(*3a)$$ is well behaved as long as $$|A|,|B|,|C| < 1$$. By continuity, $$(*3a)$$ remains valid even when $$a, b, c$$ are linear dependent.

Actually, in this special case, there is a simpler expression for $$V$$. Notice $$a,b,c$$ are determined up to a sign. Let $$\alpha = \cos^{-1}A, \beta = \cos^{-1}B, \gamma = \cos^{-1}C \quad\in (0,\pi)$$ It is always possible to choose the signs of $$b$$ and $$c$$ to make $$\alpha + \beta + \gamma = 2\pi$$. If one do that, one can simplifies $$(*3a)$$ to \bbox[border: 1px solid blue;padding: 1em]{ \begin{align} V &= \frac83 \left(\sqrt{\frac{1+A}{1-A}} + \sqrt{\frac{1+B}{1-B}} + \sqrt{\frac{1+C}{1-C}}\right)\\ &= \frac83 \left(\cot\frac{\alpha}{2} + \cot\frac{\beta}{2} + \cot\frac{\gamma}{2}\right) \end{align} }\tag{*3b} As an example, consider the case $$a, b, c$$ lies on the plane $$z = 0$$ with angle $$120^\circ$$ among them. Above formula gives $$V = \frac{8}{\sqrt{3}}$$ which can be proved with elementary means.