Consider an ellipse (say with major axis $a$ and minor axis $b$) centered at origin with a concentric circle of radius $R$. Area of the region between the circle and the ellipse is $$A = \pi R^2 – \pi a b.$$
I need to integrate a function numerically inside this region with the usual Gauss-Legendre quadrature within the classical Finite Element framework. Let the region between ellipse and the circle be divided in 4 finite elements and consider the finite element in the first quadrant. I use the following for computing the $x$, $y$ coordinates of a Gauss point in the finite element,
$$ \theta_{gauss} = \theta_0 + \eta (\theta_2 – \theta_0) $$
where, $\theta_1$ and $\theta_2$ are the angles made by the horizontal and vertical edges of the finite element considered with the $X$-axis, $\theta_0=\frac{\theta_2-\theta_1}{2}$ and $\eta$ is the first Gauss abscissa. Radius of a point that lies on the ellipse and makes an angle $\theta_{gauss}$ with the $X$-axis is simply,
$$r_e = \left(a^2 \cos^2 (\theta_{gauss}) + b^2 \sin^2 (\theta_{gauss}) \right)^{1/2}.$$
In order to place the Gauss integration points such that they are placed from the ellipse moving towards the circle, I use, $$r_0 = \frac{ R + r_e }{2} \hspace{0.1in} \mbox{and},\\
r_{gauss} = r_0 + \xi \left(\frac{R- r_e}{2}\right),$$
where, $\xi$ is the second Gauss abscissa.
Finally, the $x$ and $y$ coordinates of the Gauss integration points are now,
$$ x_{gauss} = r_{gauss} \cos (\theta_{gauss}), y_{gauss} = r_{gauss} \sin (\theta_{gauss}).$$
The problem with this though is that since I have an elliptical radius $r_e$ and a circular radius $R$ used to compute the effective radius of the Gauss point, the distribution of the integration points is not uniform. Secondly, I clearly miss a considerable area (can be seen as white space at $45^{\circ}$) and this results in erroneous results.
For this figure, a 30 x 30 Gauss rule is used in the first finite element and the red dots are the Gauss integration points (the points that are seen clustering towards the end of the ellipse do not cross over into the ellipse). Am I doing this correctly? How do I get my integration points correctly distributed in the annular region between ellipse and the circle. If this is not possible due to elliptic boundary on one side and circular on the other, is it possible to at least position the Gauss points such that the the entire annular region is accounted for (i.e. the numerically integrated area is correct).
Best Answer
I had made a mistake in the code. The integration points now appear OK and the percentage error (with 20 Gauss points in each direction) in the area is 5.9E-014. The Matlab code used for this is here in case anyone is interested.