##### https://github.com/geodynamics/citcoms
Tip revision: b5721c1
Size_does_matter.c
``````/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
* Clint Conrad, Michael Gurnis, and Eun-seo Choi.
* Copyright (C) 1994-2005, California Institute of Technology.
*
* This program is free software; you can redistribute it and/or modify
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <math.h>
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"

#ifdef ALLOW_ELLIPTICAL
double theta_g(double , struct All_variables *);
#endif

void twiddle_thumbs(struct All_variables *yawn)
//     int scratch_groin;

{ /* Do nothing, just sit back and relax.
Take it easy for a while, maybe size
doesn't matter after all. There, there
that's better. Now ... */

return; }

/*   ======================================================================
======================================================================  */

static void form_rtf_bc(int k, double x[4],
double rtf[4][9], double bc[4][4])
{
rtf[3][k] = 1.0/sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]); /* 1/r */
rtf[1][k] = acos(x[3]*rtf[3][k]); /* theta */
rtf[2][k] = myatan(x[2],x[1]); /* phi */

bc[1][1] = x[3]*cos(rtf[2][k]); /* theta */
bc[1][2] = x[3]*sin(rtf[2][k]);
bc[1][3] = -sin(rtf[1][k])/rtf[3][k];
bc[2][1] = -x[2];		/* phi basis vector */
bc[2][2] = x[1];
bc[2][3] = 0.0;
bc[3][1] = x[1]*rtf[3][k];	/*  */
bc[3][2] = x[2]*rtf[3][k];
bc[3][3] = x[3]*rtf[3][k];

return;
}

static void get_global_shape_fn_sph(struct All_variables *E,
int m, int lev, int el)
{
int i,j,k,d,e;
double jacobian;

struct Shape_function_dx LGNx;

double dxda[4][4], cof[4][4], x[4], rtf[4][9], bc[4][4];

const int dims = E->mesh.nsd;
const int ends = ENODES3D;
const int vpts = VPOINTS3D;
const int ppts = PPOINTS3D;

for(k=1;k<=vpts;k++) {       /* all of the vpoints */
for(d=1;d<=dims;d++)  {
x[d]=0.0;
for(e=1;e<=dims;e++)
dxda[d][e]=0.0;
}

for(d=1;d<=dims;d++)
for(i=1;i<=ends;i++)
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[i]]
* E->N.vpt[GNVINDEX(i,k)];

for(d=1;d<=dims;d++)
for(e=1;e<=dims;e++)
for(i=1;i<=ends;i++)
dxda[d][e] += E->X[lev][m][e][E->IEN[lev][m][el].node[i]]
* E->Nx.vpt[GNVXINDEX(d-1,i,k)];

jacobian = determinant(dxda, E->mesh.nsd);
E->GDA[lev][m][el].vpt[k] = jacobian;

for(d=1;d<=dims;d++)
for(e=1;e<=dims;e++)
cof[d][e]=cofactor(dxda,d,e,dims);

form_rtf_bc(k,x,rtf,bc);
for(j=1;j<=ends;j++)
for(d=1;d<=dims;d++)         {
LGNx.vpt[GNVXINDEX(d-1,j,k)] = 0.0;
for(e=1;e<=dims;e++)
LGNx.vpt[GNVXINDEX(d-1,j,k)] +=
E->Nx.vpt[GNVXINDEX(e-1,j,k)] *cof[e][d];

LGNx.vpt[GNVXINDEX(d-1,j,k)] /= jacobian;
}

for(j=1;j<=ends;j++)
for(d=1;d<=dims;d++)         {
E->GNX[lev][m][el].vpt[GNVXINDEX(d-1,j,k)]
= bc[d][1]*LGNx.vpt[GNVXINDEX(0,j,k)]
+ bc[d][2]*LGNx.vpt[GNVXINDEX(1,j,k)]
+ bc[d][3]*LGNx.vpt[GNVXINDEX(2,j,k)];
}
}     /* end for k */

for(k=1;k<=ppts;k++) {   /* all of the ppoints */
for(d=1;d<=dims;d++) {
x[d]=0.0;
for(e=1;e<=dims;e++)
dxda[d][e]=0.0;
}

for(d=1;d<=dims;d++)
for(i=1;i<=ends;i++)
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[i]]
* E->N.ppt[GNPINDEX(i,k)];

for(d=1;d<=dims;d++)
for(e=1;e<=dims;e++)
for(i=1;i<=ends;i++)
dxda[d][e] += E->X[lev][m][e][E->IEN[lev][m][el].node[i]]
* E->Nx.ppt[GNPXINDEX(d-1,i,k)];

jacobian = determinant(dxda,E->mesh.nsd);
E->GDA[lev][m][el].ppt[k] = jacobian;

for(d=1;d<=dims;d++)
for(e=1;e<=dims;e++)
cof[d][e]=cofactor(dxda,d,e,E->mesh.nsd);

form_rtf_bc(k,x,rtf,bc);
for(j=1;j<=ends;j++)
for(d=1;d<=dims;d++)  {
LGNx.ppt[GNPXINDEX(d-1,j,k)]=0.0;
for(e=1;e<=dims;e++)
LGNx.ppt[GNPXINDEX(d-1,j,k)] +=
E->Nx.ppt[GNPXINDEX(e-1,j,k)]*cof[e][d];
LGNx.ppt[GNPXINDEX(d-1,j,k)] /= jacobian;
}
for(j=1;j<=ends;j++)
for(d=1;d<=dims;d++)         {
E->GNX[lev][m][el].ppt[GNPXINDEX(d-1,j,k)]
= bc[d][1]*LGNx.ppt[GNPXINDEX(0,j,k)]
+ bc[d][2]*LGNx.ppt[GNPXINDEX(1,j,k)]
+ bc[d][3]*LGNx.ppt[GNPXINDEX(2,j,k)];
}

}              /* end for k int */

return;
}

void construct_shape_function_derivatives(struct All_variables *E)
{
int m, lev, el;

for (m=1; m<=E->sphere.caps_per_proc; m++)
for(lev=E->mesh.levmax; lev>=E->mesh.levmin; lev--)
for(el=1; el<=E->lmesh.NEL[lev]; el++) {
get_global_shape_fn_sph(E, m, lev, el);
}

return;
}

/*

gets r,theta,phi coordinates at the integration points

*/
void get_rtf_at_vpts(struct All_variables *E, int m, int lev, int el,
double rtf[4][9])
{
int i, k, d;
double x[4];

const int dims = E->mesh.nsd;
const int ends = ENODES3D;
const int vpts = VPOINTS3D;

for(k=1;k<=vpts;k++) {       /* all of the vpoints */
for(d=1;d<=dims;d++)
x[d]=0.0;

for(d=1;d<=dims;d++)
for(i=1;i<=ends;i++)
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[i]]
* E->N.vpt[GNVINDEX(i,k)];

rtf[3][k] = 1.0/sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]); /* 1/r */
rtf[1][k] = acos(x[3]*rtf[3][k]); /* theta */
rtf[2][k] = myatan(x[2],x[1]); /* phi */
}

return;
}

void get_rtf_at_ppts(struct All_variables *E, int m, int lev, int el,
double rtf[4][9])
{
int i, k, d;
double x[4];

const int dims = E->mesh.nsd;
const int ends = ENODES3D;
const int ppts = PPOINTS3D;

for(k=1;k<=ppts;k++) {   /* all of the ppoints */
for(d=1;d<=dims;d++)
x[d]=0.0;

for(d=1;d<=dims;d++)
for(i=1;i<=ends;i++)
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[i]]
* E->N.ppt[GNPINDEX(i,k)];

rtf[3][k] = 1.0/sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]);
rtf[1][k] = acos(x[3]*rtf[3][k]);
rtf[2][k] = myatan(x[2],x[1]);
}

return;
}

void get_side_x_cart(struct All_variables *E, double xx[4][5],
int el, int side, int m)
{
double to,fo,dxdy[4][4];
int i, node, s;
const int oned = onedvpoints[E->mesh.nsd];

#ifdef ALLOW_ELLIPTICAL
to = theta_g(E->eco[m][el].centre[1],E);
#else
to = E->eco[m][el].centre[1];
#endif

fo = E->eco[m][el].centre[2];

dxdy[1][1] = cos(to)*cos(fo);
dxdy[1][2] = cos(to)*sin(fo);
dxdy[1][3] = -sin(to);
dxdy[2][1] = -sin(fo);
dxdy[2][2] = cos(fo);
dxdy[2][3] = 0.0;
dxdy[3][1] = sin(to)*cos(fo);
dxdy[3][2] = sin(to)*sin(fo);
dxdy[3][3] = cos(to);

for(i=1;i<=oned;i++) {     /* nodes */
s = sidenodes[side][i];
node = E->ien[m][el].node[s];
xx[1][i] = E->x[m][1][node]*dxdy[1][1]
+ E->x[m][2][node]*dxdy[1][2]
+ E->x[m][3][node]*dxdy[1][3];
xx[2][i] = E->x[m][1][node]*dxdy[2][1]
+ E->x[m][2][node]*dxdy[2][2]
+ E->x[m][3][node]*dxdy[2][3];
xx[3][i] = E->x[m][1][node]*dxdy[3][1]
+ E->x[m][2][node]*dxdy[3][2]
+ E->x[m][3][node]*dxdy[3][3];
}
}

/*   ======================================================================
======================================================================  */
void construct_surf_det (struct All_variables *E)
{

int m,i,k,d,e,es,el;

double jacobian;

const int oned = onedvpoints[E->mesh.nsd];

double xx[4][5], dxda[4][4], r2;

for (m=1;m<=E->sphere.caps_per_proc;m++)
for(k=1;k<=oned;k++)    { /* all of the vpoints*/
E->surf_det[m][k] = (double *)malloc((1+E->lmesh.snel)*sizeof(double));
}

for (m=1;m<=E->sphere.caps_per_proc;m++) {
r2 = 1.0 / (E->sx[m][3][E->lmesh.elz+1] * E->sx[m][3][E->lmesh.elz+1]);

for (es=1;es<=E->lmesh.snel;es++)   {
el = es * E->lmesh.elz;
get_side_x_cart(E, xx, el, SIDE_TOP, m);

for(k=1;k<=oned;k++)    { /* all of the vpoints*/
for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e]=0.0;

for(i=1;i<=oned;i++)      /* nodes */
for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e] += xx[e][i]*E->Mx.vpt[GMVXINDEX(d-1,i,k)];

jacobian = determinant(dxda,E->mesh.nsd-1);

/* scale the jacobian so that it is defined on a unit sphere */
E->surf_det[m][k][es] = jacobian * r2;
}
}
}
return;
}

/*   ======================================================================
surface (6 sides) determinant of boundary element
======================================================================  */
void construct_bdry_det(struct All_variables *E)
{

int m,i,k,d,e,es,el,side;

double jacobian;

const int oned = onedvpoints[E->mesh.nsd];

double xx[4][5],dxda[4][4];

for (m=1;m<=E->sphere.caps_per_proc;m++)
for (side=SIDE_BEGIN; side<=SIDE_END; side++)
for(d=1; d<=oned; d++)
E->boundary.det[m][side][d] = (double *)malloc((1+E->boundary.nel)*sizeof(double));

for (m=1;m<=E->sphere.caps_per_proc;m++)
for (es=1;es<=E->boundary.nel;es++) {
el = E->boundary.element[m][es];

for (side=SIDE_BEGIN; side<=SIDE_END; side++) {
get_side_x_cart(E, xx, el, side, m);

for(k=1;k<=oned;k++) { /* all of the vpoints*/

for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e]=0.0;

for(i=1;i<=oned;i++) /* nodes */
for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e] += xx[sidedim[side][e]][i]*E->Mx.vpt[GMVXINDEX(d-1,i,k)];

jacobian = determinant(dxda,E->mesh.nsd-1);
E->boundary.det[m][side][k][es] = jacobian;
}

/*
fprintf(stderr, "Boundary det: %d %d- %e %e %e %e; sum = %e\n", el, side,
E->boundary.det[m][side][1][es],
E->boundary.det[m][side][2][es],
E->boundary.det[m][side][3][es],
E->boundary.det[m][side][4][es],
E->boundary.det[m][side][1][es]+
E->boundary.det[m][side][2][es]+
E->boundary.det[m][side][3][es]+
E->boundary.det[m][side][4][es]);
*/
}

}
}

/*   ======================================================================
======================================================================  */
void get_global_1d_shape_fn(struct All_variables *E, int el, struct Shape_function1 *GM, struct Shape_function1_dA *dGammax, int top, int m)
{
int ii,i,k,d,e;

double jacobian;

const int oned = onedvpoints[E->mesh.nsd];

double xx[4][5],dxda[4][4];

for (ii=0;ii<=top;ii++)   {   /* ii=0 for bottom and ii=1 for top */

get_side_x_cart(E, xx, el, ii+1, m);

for(k=1;k<=oned;k++)    { /* all of the vpoints*/
for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e]=0.0;

for(i=1;i<=oned;i++)      /* nodes */
for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e] += xx[e][i]*E->Mx.vpt[GMVXINDEX(d-1,i,k)];

jacobian = determinant(dxda,E->mesh.nsd-1);
dGammax->vpt[GMVGAMMA(ii,k)] = jacobian;
}
}

return;
}

/*   ======================================================================
======================================================================  */
void get_global_1d_shape_fn_L(struct All_variables *E, int el, struct Shape_function1 *GM, struct Shape_function1_dA *dGammax, int top, int m)
{
int ii,i,k,d,e,node;

double jacobian;

const int oned = onedvpoints[E->mesh.nsd];

double to,fo,xx[4][5],dxdy[4][4],dxda[4][4],cof[4][4];

#ifdef ALLOW_ELLIPTICAL
to = theta_g(E->eco[m][el].centre[1],E);
#else
to = E->eco[m][el].centre[1];
#endif
fo = E->eco[m][el].centre[2];

dxdy[1][1] = cos(to)*cos(fo);
dxdy[1][2] = cos(to)*sin(fo);
dxdy[1][3] = -sin(to);
dxdy[2][1] = -sin(fo);
dxdy[2][2] = cos(fo);
dxdy[2][3] = 0.0;
dxdy[3][1] = sin(to)*cos(fo);
dxdy[3][2] = sin(to)*sin(fo);
dxdy[3][3] = cos(to);

for (ii=0;ii<=top;ii++)   {   /* ii=0 for bottom and ii=1 for top */

for(i=1;i<=oned;i++) {     /* nodes */
e = i+ii*oned;
node = E->ien[m][el].node[e];
xx[1][i] = E->x[m][1][node]*dxdy[1][1]
+ E->x[m][2][node]*dxdy[1][2]
+ E->x[m][3][node]*dxdy[1][3];
xx[2][i] = E->x[m][1][node]*dxdy[2][1]
+ E->x[m][2][node]*dxdy[2][2]
+ E->x[m][3][node]*dxdy[2][3];
xx[3][i] = E->x[m][1][node]*dxdy[3][1]
+ E->x[m][2][node]*dxdy[3][2]
+ E->x[m][3][node]*dxdy[3][3];
}

for(k=1;k<=oned;k++)    { /* all of the vpoints*/

for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e]=0.0;

for(i=1;i<=oned;i++)      /* nodes */
for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e] += xx[e][i]*E->Lx.vpt[GMVXINDEX(d-1,i,k)];

jacobian = determinant(dxda,E->mesh.nsd-1);
dGammax->vpt[GMVGAMMA(ii,k)] = jacobian;
}
}

return;
}

/*   ======================================================================
For calculating pressure boundary term --- Choi, 11/13/02
======================================================================  */
void get_global_side_1d_shape_fn(struct All_variables *E, int el, struct Shape_function1 *GM, struct Shape_function1_dx *GMx, struct Shape_function_side_dA *dGamma, int side, int m)
{
int i,k,d,e;

double jacobian;

const int oned = onedvpoints[E->mesh.nsd];
double xx[4][5],dxda[4][4];

get_side_x_cart(E, xx, el, side, m);

for(k=1;k<=oned;k++)    {

for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++)
dxda[d][e]=0.0;

for(i=1;i<=oned;i++) {
for(d=1;d<=E->mesh.nsd-1;d++)
for(e=1;e<=E->mesh.nsd-1;e++) {
dxda[d][e] += xx[sidedim[side][e]][i]*E->Mx.vpt[GMVXINDEX(d-1,i,k)];
}
}

jacobian = determinant(dxda,E->mesh.nsd-1);
dGamma->vpt[k] = jacobian;
}

return;
}

/* ====================================================   */

void construct_c3x3matrix_el (struct All_variables *E,int el,struct CC *cc,
struct CCX *ccx,int lev,int m,int pressure)
{
int a,i,j,k,d,lnode;
double x[4],u[4][4],ux[3][4][4],ua[4][4];
double costt,cosff,sintt,sinff,rr,tt,ff;

const int dims=E->mesh.nsd;
const int ends=enodes[dims];
const int vpts=vpoints[dims];
const int ppts=ppoints[dims];

if (pressure==0)           {
for(k=1;k<=vpts;k++)           {       /* all of the vpoints */
for(d=1;d<=dims;d++)
x[d]=0.0;

for(d=1;d<=dims;d++)
for(a=1;a<=ends;a++)
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[a]]
*E->N.vpt[GNVINDEX(a,k)];

rr = sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]);

ff = myatan(x[2],x[1]);
#ifdef ALLOW_ELLIPTICAL
tt = theta_g(acos(x[3]/rr),E);
#else
tt = acos(x[3]/rr);
#endif
costt = cos(tt);
cosff = cos(ff);
sintt = sin(tt);
sinff = sin(ff);

u[1][1] = costt*cosff; u[1][2] = costt*sinff;  u[1][3] =-sintt;
u[2][1] =-sinff;       u[2][2] = cosff;        u[2][3] = 0.0;
u[3][1] = sintt*cosff; u[3][2] = sintt*sinff;  u[3][3] = costt;

ux[1][1][1] =-sintt*cosff;  ux[1][1][2] =-sintt*sinff;  ux[1][1][3] =-costt;
ux[2][1][1] =-costt*sinff;  ux[2][1][2] = costt*cosff;  ux[2][1][3] =0.0;
ux[1][2][1] =0.0;           ux[1][2][2] = 0.0;          ux[1][2][3] =0.0;
ux[2][2][1] =-cosff;        ux[2][2][2] =-sinff;        ux[2][2][3] =0.0;
ux[1][3][1] = costt*cosff;  ux[1][3][2] = costt*sinff;  ux[1][3][3] =-sintt;
ux[2][3][1] =-sintt*sinff;  ux[2][3][2] = sintt*cosff;  ux[2][3][3] =0.0;

for(a=1;a<=ends;a++)   {
tt = E->SX[lev][m][1][E->IEN[lev][m][el].node[a]];
ff = E->SX[lev][m][2][E->IEN[lev][m][el].node[a]];
costt = cos(tt);
cosff = cos(ff);
sintt = sin(tt);
sinff = sin(ff);

ua[1][1] = costt*cosff; ua[1][2] = costt*sinff;  ua[1][3] =-sintt;
ua[2][1] =-sinff;       ua[2][2] = cosff;        ua[2][3] = 0.0;
ua[3][1] = sintt*cosff; ua[3][2] = sintt*sinff;  ua[3][3] = costt;

for (i=1;i<=dims;i++)
for (j=1;j<=dims;j++)   {
cc->vpt[BVINDEX(i,j,a,k)] =
ua[j][1]*u[i][1]+ua[j][2]*u[i][2]+ua[j][3]*u[i][3];
ccx->vpt[BVXINDEX(i,j,1,a,k)] =
ua[j][1]*ux[1][i][1]+ua[j][2]*ux[1][i][2]+ua[j][3]*ux[1][i][3];
ccx->vpt[BVXINDEX(i,j,2,a,k)] =
ua[j][1]*ux[2][i][1]+ua[j][2]*ux[2][i][2]+ua[j][3]*ux[2][i][3];
}
}      /* end for local node */

}        /* end for int points */
}        /* end if */

else if (pressure)  {

for(k=1;k<=ppts;k++)           {       /* all of the ppoints */
for(d=1;d<=dims;d++)
x[d]=0.0;

for(d=1;d<=dims;d++)
for(a=1;a<=ends;a++)
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[a]]
*E->N.ppt[GNPINDEX(a,k)];

rr = sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]);

ff = myatan(x[2],x[1]);
#ifdef ALLOW_ELLIPTICAL
tt = theta_g(acos(x[3]/rr),E);
#else
tt = acos(x[3]/rr);
#endif
costt = cos(tt);
cosff = cos(ff);
sintt = sin(tt);
sinff = sin(ff);

u[1][1] = costt*cosff; u[1][2] = costt*sinff;  u[1][3] =-sintt;
u[2][1] =-sinff;       u[2][2] = cosff;        u[2][3] = 0.0;
u[3][1] = sintt*cosff; u[3][2] = sintt*sinff;  u[3][3] = costt;

ux[1][1][1] =-sintt*cosff;  ux[1][1][2] =-sintt*sinff;  ux[1][1][3] =-costt;
ux[2][1][1] =-costt*sinff;  ux[2][1][2] = costt*cosff;  ux[2][1][3] =0.0;
ux[1][2][1] =0.0;           ux[1][2][2] = 0.0;          ux[1][2][3] =0.0;
ux[2][2][1] =-cosff;        ux[2][2][2] =-sinff;        ux[2][2][3] =0.0;
ux[1][3][1] = costt*cosff;  ux[1][3][2] = costt*sinff;  ux[1][3][3] =-sintt;
ux[2][3][1] =-sintt*sinff;  ux[2][3][2] = sintt*cosff;  ux[2][3][3] =0.0;

for(a=1;a<=ends;a++)   {

lnode = E->IEN[lev][m][el].node[a];
sintt = E->SinCos[lev][m][0][lnode];
sinff = E->SinCos[lev][m][1][lnode];
costt = E->SinCos[lev][m][2][lnode];
cosff = E->SinCos[lev][m][3][lnode];

ua[1][1] = costt*cosff; ua[1][2] = costt*sinff;  ua[1][3] =-sintt;
ua[2][1] =-sinff;       ua[2][2] = cosff;        ua[2][3] = 0.0;
ua[3][1] = sintt*cosff; ua[3][2] = sintt*sinff;  ua[3][3] = costt;

for (i=1;i<=dims;i++)
for (j=1;j<=dims;j++)   {
cc->ppt[BPINDEX(i,j,a,k)] =
ua[j][1]*u[i][1]+ua[j][2]*u[i][2]+ua[j][3]*u[i][3];
ccx->ppt[BPXINDEX(i,j,1,a,k)] =
ua[j][1]*ux[1][i][1]+ua[j][2]*ux[1][i][2]+ua[j][3]*ux[1][i][3];
ccx->ppt[BPXINDEX(i,j,2,a,k)] =
ua[j][1]*ux[2][i][1]+ua[j][2]*ux[2][i][2]+ua[j][3]*ux[2][i][3];
}

}      /* end for local node */

}        /* end for int points */

}         /* end if pressure  */

return;
}

void construct_side_c3x3matrix_el(struct All_variables *E,int el,
struct CC *cc,struct CCX *ccx,
int lev,int m,int pressure,int side)
{
int a,aa,i,j,k,d,lnode;
double x[4],u[4][4],ux[3][4][4],ua[4][4];
double costt,cosff,sintt,sinff,rr,tt,ff;

const int dims=E->mesh.nsd;
const int ends=enodes[dims-1];
const int vpts=onedvpoints[dims];
const int ppts=ppoints[dims];

if(pressure==0) {
for(k=1;k<=vpts;k++) {       /* all of the vpoints */
for(d=1;d<=dims;d++)
x[d]=0.0;
for(d=1;d<=dims;d++)
for(aa=1;aa<=ends;aa++) {
a=sidenodes[side][aa];
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[a]]
*E->M.vpt[GMVINDEX(aa,k)];

}

rr = sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]);
tt = acos(x[3]/rr);
ff = myatan(x[2],x[1]);
#ifdef ALLOW_ELLIPTICAL
tt = theta_g(acos(x[3]/rr),E);
#else
tt = acos(x[3]/rr);
#endif

costt = cos(tt);
cosff = cos(ff);
sintt = sin(tt);
sinff = sin(ff);

u[1][1] = costt*cosff; u[1][2] = costt*sinff;  u[1][3] =-sintt;
u[2][1] =-sinff;       u[2][2] = cosff;        u[2][3] = 0.0;
u[3][1] = sintt*cosff; u[3][2] = sintt*sinff;  u[3][3] = costt;

ux[1][1][1] =-sintt*cosff;  ux[1][1][2] =-sintt*sinff;  ux[1][1][3] =-costt;
ux[2][1][1] =-costt*sinff;  ux[2][1][2] = costt*cosff;  ux[2][1][3] =0.0;
ux[1][2][1] =0.0;           ux[1][2][2] = 0.0;          ux[1][2][3] =0.0;
ux[2][2][1] =-cosff;        ux[2][2][2] =-sinff;        ux[2][2][3] =0.0;
ux[1][3][1] = costt*cosff;  ux[1][3][2] = costt*sinff;  ux[1][3][3] =-sintt;
ux[2][3][1] =-sintt*sinff;  ux[2][3][2] = sintt*cosff;  ux[2][3][3] =0.0;

for(aa=1;aa<=ends;aa++) {
a=sidenodes[side][aa];

lnode = E->IEN[lev][m][el].node[a];
sintt = E->SinCos[lev][m][0][lnode];
sinff = E->SinCos[lev][m][1][lnode];
costt = E->SinCos[lev][m][2][lnode];
cosff = E->SinCos[lev][m][3][lnode];

ua[1][1] = costt*cosff; ua[1][2] = costt*sinff;  ua[1][3] =-sintt;
ua[2][1] =-sinff;       ua[2][2] = cosff;        ua[2][3] = 0.0;
ua[3][1] = sintt*cosff; ua[3][2] = sintt*sinff;  ua[3][3] = costt;

for (i=1;i<=dims;i++)
for (j=1;j<=dims;j++)   {
cc->vpt[BVINDEX(i,j,a,k)] =
ua[j][1]*u[i][1]+ua[j][2]*u[i][2]+ua[j][3]*u[i][3];
ccx->vpt[BVXINDEX(i,j,1,a,k)] =
ua[j][1]*ux[1][i][1]+ua[j][2]*ux[1][i][2]+ua[j][3]*ux[1][i][3];
ccx->vpt[BVXINDEX(i,j,2,a,k)] =
ua[j][1]*ux[2][i][1]+ua[j][2]*ux[2][i][2]+ua[j][3]*ux[2][i][3];
}

}      /* end for local node */
}        /* end for int points */
}    /* end if */
else {
for(k=1;k<=ppts;k++) {       /* all of the ppoints */
for(d=1;d<=E->mesh.nsd;d++)
x[d]=0.0;
for(a=1;a<=ends;a++) {
aa=sidenodes[side][a];
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[aa]]
*E->M.ppt[GMPINDEX(a,k)];
}
rr = sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]);
ff = myatan(x[2],x[1]);

#ifdef ALLOW_ELLIPTICAL
tt = theta_g(acos(x[3]/rr),E);
#else
tt = acos(x[3]/rr);
#endif
costt = cos(tt);
cosff = cos(ff);
sintt = sin(tt);
sinff = sin(ff);

u[1][1] = costt*cosff; u[1][2] = costt*sinff;  u[1][3] =-sintt;
u[2][1] =-sinff;       u[2][2] = cosff;        u[2][3] = 0.0;
u[3][1] = sintt*cosff; u[3][2] = sintt*sinff;  u[3][3] = costt;

ux[1][1][1] =-sintt*cosff;  ux[1][1][2] =-sintt*sinff;  ux[1][1][3] =-costt;
ux[2][1][1] =-costt*sinff;  ux[2][1][2] = costt*cosff;  ux[2][1][3] =0.0;
ux[1][2][1] =0.0;           ux[1][2][2] = 0.0;          ux[1][2][3] =0.0;
ux[2][2][1] =-cosff;        ux[2][2][2] =-sinff;        ux[2][2][3] =0.0;
ux[1][3][1] = costt*cosff;  ux[1][3][2] = costt*sinff;  ux[1][3][3] =-sintt;
ux[2][3][1] =-sintt*sinff;  ux[2][3][2] = sintt*cosff;  ux[2][3][3] =0.0;

for(a=1;a<=ends;a++)   {
aa=sidenodes[side][a];

lnode = E->IEN[lev][m][el].node[aa];
sintt = E->SinCos[lev][m][0][lnode];
sinff = E->SinCos[lev][m][1][lnode];
costt = E->SinCos[lev][m][2][lnode];
cosff = E->SinCos[lev][m][3][lnode];

ua[1][1] = costt*cosff; ua[1][2] = costt*sinff;  ua[1][3] =-sintt;
ua[2][1] =-sinff;       ua[2][2] = cosff;        ua[2][3] = 0.0;
ua[3][1] = sintt*cosff; ua[3][2] = sintt*sinff;  ua[3][3] = costt;

for (i=1;i<=E->mesh.nsd;i++) {
for (j=1;j<=E->mesh.nsd;j++) {
cc->ppt[BPINDEX(i,j,a,k)] =
ua[j][1]*u[i][1]+ua[j][2]*u[i][2]+ua[j][3]*u[i][3];
ccx->ppt[BPXINDEX(i,j,1,a,k)] =
ua[j][1]*ux[1][i][1]+ua[j][2]*ux[1][i][2]+ua[j][3]*ux[1][i][3];
ccx->ppt[BPXINDEX(i,j,2,a,k)] =
ua[j][1]*ux[2][i][1]+ua[j][2]*ux[2][i][2]+ua[j][3]*ux[2][i][3];
}
}
}      /* end for local node */
}      /* end for int points */
}      /* end if pressure  */

return;
}

/* ======================================= */
void construct_c3x3matrix(struct All_variables *E)
{
int m,a,i,j,k,d,es,el,nel_surface,lev,lnode;
double x[4],u[4][4],ux[3][4][4],ua[4][4];
double costt,cosff,sintt,sinff,rr,tt,ff;

const int dims=E->mesh.nsd;
const int ends=enodes[dims];
const int vpts=vpoints[dims];
const int ppts=ppoints[dims];

for (lev=E->mesh.gridmin;lev<=E->mesh.gridmax;lev++)
for (m=1;m<=E->sphere.caps_per_proc;m++)   {
nel_surface = E->lmesh.NEL[lev]/E->lmesh.ELZ[lev];
for (es=1;es<=nel_surface;es++)        {

el = es*E->lmesh.ELZ[lev];

for(k=1;k<=vpts;k++)           {       /* all of the vpoints */
for(d=1;d<=dims;d++)
x[d]=0.0;

for(d=1;d<=dims;d++)
for(a=1;a<=ends;a++)
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[a]]
*E->N.vpt[GNVINDEX(a,k)];

rr = sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]);
ff = myatan(x[2],x[1]);
#ifdef ALLOW_ELLIPTICAL
tt = theta_g(acos(x[3]/rr),E);
#else
tt = acos(x[3]/rr);
#endif

costt = cos(tt);
cosff = cos(ff);
sintt = sin(tt);
sinff = sin(ff);

u[1][1] = costt*cosff; u[1][2] = costt*sinff;  u[1][3] =-sintt;
u[2][1] =-sinff;       u[2][2] = cosff;        u[2][3] = 0.0;
u[3][1] = sintt*cosff; u[3][2] = sintt*sinff;  u[3][3] = costt;

ux[1][1][1] =-sintt*cosff;  ux[1][1][2] =-sintt*sinff;  ux[1][1][3] =-costt;
ux[2][1][1] =-costt*sinff;  ux[2][1][2] = costt*cosff;  ux[2][1][3] =0.0;
ux[1][2][1] =0.0;           ux[1][2][2] = 0.0;          ux[1][2][3] =0.0;
ux[2][2][1] =-cosff;        ux[2][2][2] =-sinff;        ux[2][2][3] =0.0;
ux[1][3][1] = costt*cosff;  ux[1][3][2] = costt*sinff;  ux[1][3][3] =-sintt;
ux[2][3][1] =-sintt*sinff;  ux[2][3][2] = sintt*cosff;  ux[2][3][3] =0.0;

for(a=1;a<=ends;a++)   {

lnode = E->IEN[lev][m][el].node[a];
sintt = E->SinCos[lev][m][0][lnode];
sinff = E->SinCos[lev][m][1][lnode];
costt = E->SinCos[lev][m][2][lnode];
cosff = E->SinCos[lev][m][3][lnode];

ua[1][1] = costt*cosff; ua[1][2] = costt*sinff;  ua[1][3] =-sintt;
ua[2][1] =-sinff;       ua[2][2] = cosff;        ua[2][3] = 0.0;
ua[3][1] = sintt*cosff; ua[3][2] = sintt*sinff;  ua[3][3] = costt;

for (i=1;i<=dims;i++)
for (j=1;j<=dims;j++)   {
E->CC[lev][m][es].vpt[BVINDEX(i,j,a,k)] =
ua[j][1]*u[i][1]+ua[j][2]*u[i][2]+ua[j][3]*u[i][3];
E->CCX[lev][m][es].vpt[BVXINDEX(i,j,1,a,k)] =
ua[j][1]*ux[1][i][1]+ua[j][2]*ux[1][i][2]+ua[j][3]*ux[1][i][3];
E->CCX[lev][m][es].vpt[BVXINDEX(i,j,2,a,k)] =
ua[j][1]*ux[2][i][1]+ua[j][2]*ux[2][i][2]+ua[j][3]*ux[2][i][3];
}
}      /* end for local node */

}        /* end for int points */

for(k=1;k<=ppts;k++)           {       /* all of the ppoints */
for(d=1;d<=dims;d++)
x[d]=0.0;

for(d=1;d<=dims;d++)
for(a=1;a<=ends;a++)
x[d] += E->X[lev][m][d][E->IEN[lev][m][el].node[a]]
*E->N.ppt[GNPINDEX(a,k)];

rr = sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]);
ff = myatan(x[2],x[1]);
#ifdef ALLOW_ELLIPTICAL
tt = theta_g(acos(x[3]/rr),E);
#else
tt = acos(x[3]/rr);
#endif

costt = cos(tt);
cosff = cos(ff);
sintt = sin(tt);
sinff = sin(ff);

u[1][1] = costt*cosff; u[1][2] = costt*sinff;  u[1][3] =-sintt;
u[2][1] =-sinff;       u[2][2] = cosff;        u[2][3] = 0.0;
u[3][1] = sintt*cosff; u[3][2] = sintt*sinff;  u[3][3] = costt;

ux[1][1][1] =-sintt*cosff;  ux[1][1][2] =-sintt*sinff;  ux[1][1][3] =-costt;
ux[2][1][1] =-costt*sinff;  ux[2][1][2] = costt*cosff;  ux[2][1][3] =0.0;
ux[1][2][1] =0.0;           ux[1][2][2] = 0.0;          ux[1][2][3] =0.0;
ux[2][2][1] =-cosff;        ux[2][2][2] =-sinff;        ux[2][2][3] =0.0;
ux[1][3][1] = costt*cosff;  ux[1][3][2] = costt*sinff;  ux[1][3][3] =-sintt;
ux[2][3][1] =-sintt*sinff;  ux[2][3][2] = sintt*cosff;  ux[2][3][3] =0.0;

for(a=1;a<=ends;a++)   {

lnode = E->IEN[lev][m][el].node[a];
sintt = E->SinCos[lev][m][0][lnode];
sinff = E->SinCos[lev][m][1][lnode];
costt = E->SinCos[lev][m][2][lnode];
cosff = E->SinCos[lev][m][3][lnode];

ua[1][1] = costt*cosff; ua[1][2] = costt*sinff;  ua[1][3] =-sintt;
ua[2][1] =-sinff;       ua[2][2] = cosff;        ua[2][3] = 0.0;
ua[3][1] = sintt*cosff; ua[3][2] = sintt*sinff;  ua[3][3] = costt;

for (i=1;i<=dims;i++)
for (j=1;j<=dims;j++)   {
E->CC[lev][m][es].ppt[BPINDEX(i,j,a,k)] =
ua[j][1]*u[i][1]+ua[j][2]*u[i][2]+ua[j][3]*u[i][3];
E->CCX[lev][m][es].ppt[BPXINDEX(i,j,1,a,k)] =
ua[j][1]*ux[1][i][1]+ua[j][2]*ux[1][i][2]+ua[j][3]*ux[1][i][3];
E->CCX[lev][m][es].ppt[BPXINDEX(i,j,2,a,k)] =
ua[j][1]*ux[2][i][1]+ua[j][2]*ux[2][i][2]+ua[j][3]*ux[2][i][3];
}

}      /* end for local node */

}        /* end for int points */

}         /* end for es */
}           /* end for m */

return;
}

/*  ==========================================
construct the lumped mass matrix. The full
matrix is the FE integration of the density
field. The lumped version is the diagonal
matrix obtained by letting the shape function
Na be delta(a,b)
========================================== */

void mass_matrix(struct All_variables *E)
{
int m,node,i,nint,e,lev;
int n[9], nz;
double area,centre[4],temp[9],temp2[9],dx1,dx2,dx3;

const int vpts=vpoints[E->mesh.nsd];

/* ECO .size can also be defined here */

for(lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++)  {
for (m=1;m<=E->sphere.caps_per_proc;m++)   {

for(node=1;node<=E->lmesh.NNO[lev];node++)
E->MASS[lev][m][node] = 0.0;

for(e=1;e<=E->lmesh.NEL[lev];e++)  {

area = centre[1] = centre[2] = centre[3] = 0.0;

for(node=1;node<=enodes[E->mesh.nsd];node++)
n[node] = E->IEN[lev][m][e].node[node];

for(i=1;i<=E->mesh.nsd;i++)  {
for(node=1;node<=enodes[E->mesh.nsd];node++)
centre[i] += E->X[lev][m][i][n[node]];

centre[i] = centre[i]/enodes[E->mesh.nsd];
}     /* end for i */

dx3 = sqrt(centre[1]*centre[1]+centre[2]*centre[2]+centre[3]*centre[3]);
dx1 = acos( centre[3]/dx3 );
dx2 = myatan(centre[2],centre[1]);

/* center of this element in the spherical coordinate */
E->ECO[lev][m][e].centre[1] = dx1;
E->ECO[lev][m][e].centre[2] = dx2;
E->ECO[lev][m][e].centre[3] = dx3;

/* delta(theta) of this element */
dx1 = citmax( fabs(E->SX[lev][m][1][n[3]]-E->SX[lev][m][1][n[1]]),
fabs(E->SX[lev][m][1][n[2]]-E->SX[lev][m][1][n[4]]) );

/* length of this element in the theta-direction */
E->ECO[lev][m][e].size[1] = dx1*E->ECO[lev][m][e].centre[3];

/* delta(phi) of this element */
dx1 = fabs(E->SX[lev][m][2][n[3]]-E->SX[lev][m][2][n[1]]);
if (dx1>M_PI)
dx1 = citmin(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) + 2.0*M_PI -
citmax(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) ;

dx2 = fabs(E->SX[lev][m][2][n[2]]-E->SX[lev][m][2][n[4]]);
if (dx2>M_PI)
dx2 = citmin(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) + 2.0*M_PI -
citmax(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) ;

dx2 = citmax(dx1,dx2);

/* length of this element in the phi-direction */
E->ECO[lev][m][e].size[2] = dx2*E->ECO[lev][m][e].centre[3]
*sin(E->ECO[lev][m][e].centre[1]);

/* delta(radius) of this element */
dx3 = 0.25*(fabs(E->SX[lev][m][3][n[5]]+E->SX[lev][m][3][n[6]]
+E->SX[lev][m][3][n[7]]+E->SX[lev][m][3][n[8]]
-E->SX[lev][m][3][n[1]]-E->SX[lev][m][3][n[2]]
-E->SX[lev][m][3][n[3]]-E->SX[lev][m][3][n[4]]));

/* length of this element in the radius-direction */
E->ECO[lev][m][e].size[3] = dx3;

/* volume (area in 2D) of this element */
for(nint=1;nint<=vpts;nint++)
area += g_point[nint].weight[E->mesh.nsd-1] * E->GDA[lev][m][e].vpt[nint];
E->ECO[lev][m][e].area = area;

for(node=1;node<=enodes[E->mesh.nsd];node++)  {
temp[node] = 0.0;
for(nint=1;nint<=vpts;nint++)
temp[node] += E->GDA[lev][m][e].vpt[nint]*g_point[nint].weight[E->mesh.nsd-1]
*E->N.vpt[GNVINDEX(node,nint)];       /* int Na dV */
}

for(node=1;node<=enodes[E->mesh.nsd];node++)
E->MASS[lev][m][E->IEN[lev][m][e].node[node]] += temp[node];

/* weight of each node, equivalent to pmass in ConMan */
for(node=1;node<=enodes[E->mesh.nsd];node++)
E->TWW[lev][m][e].node[node] = temp[node];

} /* end of ele*/

} /* end of for m */

if(lev == E->mesh.levmax)
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.NNO[lev];node++)
E->NMass[m][node] = E->MASS[lev][m][node];

if (E->control.NMULTIGRID||E->mesh.levmax==lev)
(E->exchange_node_d)(E,E->MASS[lev],lev);

for (m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.NNO[lev];node++)
E->MASS[lev][m][node] = 1.0/E->MASS[lev][m][node];

} /* end of for lev */

for (m=1;m<=E->sphere.caps_per_proc;m++) {

for(node=1;node<=E->lmesh.nno;node++)
E->TMass[m][node] = 0.0;

for(e=1;e<=E->lmesh.nel;e++)  {
for(node=1;node<=enodes[E->mesh.nsd];node++) {
temp[node] = 0.0;
nz = ((E->ien[m][e].node[node]-1) % E->lmesh.noz) + 1;
for(nint=1;nint<=vpts;nint++)
temp[node] += E->refstate.rho[nz]
* E->refstate.heat_capacity[nz]
* E->gDA[m][e].vpt[nint]
* g_point[nint].weight[E->mesh.nsd-1]
* E->N.vpt[GNVINDEX(node,nint)];
}

/* lumped mass matrix, equivalent to tmass in ConMan */
for(node=1;node<=enodes[E->mesh.nsd];node++)
E->TMass[m][E->ien[m][e].node[node]] += temp[node];

} /* end of for e */
} /* end of for m */

(E->exchange_node_d)(E,E->TMass,E->mesh.levmax);
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.nno;node++)
E->TMass[m][node] = 1.0 / E->TMass[m][node];

/* compute volume of this processor mesh and the whole mesh */
E->lmesh.volume = 0;
E->mesh.volume = 0;

for (m=1;m<=E->sphere.caps_per_proc;m++)
for(e=1;e<=E->lmesh.nel;e++)
E->lmesh.volume += E->eco[m][e].area;

MPI_Allreduce(&E->lmesh.volume, &E->mesh.volume, 1, MPI_DOUBLE,
MPI_SUM, E->parallel.world);

if (E->control.verbose)  {
fprintf(E->fp_out, "rank=%d my_volume=%e total_volume=%e\n",
E->parallel.me, E->lmesh.volume, E->mesh.volume);

for(lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++)  {
fprintf(E->fp_out,"output_mass lev=%d\n",lev);
for (m=1;m<=E->sphere.caps_per_proc;m++)   {
fprintf(E->fp_out,"m=%d %d \n",E->sphere.capid[m],m);
for(e=1;e<=E->lmesh.NEL[lev];e++)
fprintf(E->fp_out,"%d %g \n",e,E->ECO[lev][m][e].area);
for (node=1;node<=E->lmesh.NNO[lev];node++)
fprintf(E->fp_out,"Mass[%d]= %g \n",node,E->MASS[lev][m][node]);
}
}

for (m=1;m<=E->sphere.caps_per_proc;m++)   {
fprintf(E->fp_out,"m=%d %d \n",E->sphere.capid[m],m);
for (node=1;node<=E->lmesh.nno;node++)
fprintf(E->fp_out,"TMass[%d]= %g \n",node,E->TMass[m][node]);
}
fflush(E->fp_out);
}

return;
}

/* version */
/* \$Id\$ */

/* End of file  */
``````