[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]

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=...)
        {
           Point<dim> x = fe_values.quadrature_point(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/