# [aspect-devel] How to grab the radius of every interior sphere made up from cell faces in the domain and some MPI questions

Wolfgang Bangerth bangerth at tamu.edu
Sun Oct 25 09:01:47 PDT 2015

Shangxin & all,
I am also a bit late to the discussion, but let me give me \$0.02 anyway:

I think you started this by finding an algorithm that is well suited to one
kind of numerical scheme (namely uniform meshes in spherical coordinates) and
then trying to fit that into a code that implements a different method. This
sometimes works, but it often leads to algorithms that are unnecessarily
complicated and inefficient.

Instead, let's start from what finite element codes do really well on general
unstructured meshes (such as the ones used by ASPECT): computing volume
integrals. In your case, if you subdivide the Earth into layers, then
expanding the density into spherical harmonics requires you to compute
integrals of the form

r_{lmi} =
\int_{r_i}^{r_{i+1}} \int \int  drho(r,phi,theta) Y_{lm} dtheta dphi dr

which you can quite easily implement as follows:
{
for (cell=...)
if (cell->is_locally_owned())
for (q_point=...)
{
double r = x.norm();
int i = ...find index of layer for radius r...;

rho = ...get rho at x;
drho = rho - rho_reference(r);

for (l=...)
for (m=...)
r[l][m][i] += drho * Y(l,m,x) * fe_values.JxW(q_point);

Y = Utilities::MPI::sum(Y, mpi_communicator);
}

I believe this is in essence what Ian is describing.

Now the question is how you get your layer subdivision r_i ... r_{i+1}. You
tried to derive it from the mesh itself, but that's a brittle approach that,
as you noticed, is not well suites to parallel computing or adaptive meshes.
My suggestion would be to simply specify a subdivision yourself, either
equidistant or otherwise, with a number of layers provided by the user.

I hope this makes sense! Best
Wolfgang

--
------------------------------------------------------------------------
Wolfgang Bangerth               email:            bangerth at math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/