https://github.com/geodynamics/citcoms
Revision df821c06952128812763ec5b7a944cd178ca5ed3 authored by Thorsten Becker on 08 December 2010, 18:43:58 UTC, committed by Thorsten Becker on 08 December 2010, 18:43:58 UTC
GMT/ggrd compile before).
1 parent cf7495e
Tip revision: df821c06952128812763ec5b7a944cd178ca5ed3 authored by Thorsten Becker on 08 December 2010, 18:43:58 UTC
Added definition of TRUE (1) and FALSE (0) in case undefined (was defined for
Added definition of TRUE (1) and FALSE (0) in case undefined (was defined for
Tip revision: df821c0
Solver_multigrid.c
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<LicenseText>
*
* 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
* it under the terms of the GNU General Public License as published by
* 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
*
*</LicenseText>
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include "element_definitions.h"
#include "global_defs.h"
#include <math.h>
void set_mg_defaults(E)
struct All_variables *E;
{ void assemble_forces_iterative();
void solve_constrained_flow_iterative();
void mg_allocate_vars();
E->control.NMULTIGRID = 1;
E->build_forcing_term = assemble_forces_iterative;
E->solve_stokes_problem = solve_constrained_flow_iterative;
E->solver_allocate_vars = mg_allocate_vars;
return;
}
void mg_allocate_vars(E)
struct All_variables *E;
{
return;
}
/* =====================================================
Function to inject data from fine to coarse grid (i.e.
just dropping values at shared grid points.
===================================================== */
void inject_scalar(E,start_lev,AU,AD)
struct All_variables *E;
int start_lev;
float **AU,**AD; /* data on upper/lower mesh */
{
int i,m,el,node_coarse,node_fine,sl_minus,eqn,eqn_coarse;
void gather_to_1layer_node ();
const int dims = E->mesh.nsd;
const int ends = enodes[dims];
if(start_lev == E->mesh.levmin) {
fprintf(E->fp,"Warning, attempting to project below lowest level\n");
return;
}
sl_minus = start_lev-1;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=E->lmesh.NEL[sl_minus];el++)
for(i=1;i<=ends;i++) {
node_coarse = E->IEN[sl_minus][m][el].node[i];
node_fine=E->IEN[start_lev][m][E->EL[sl_minus][m][el].sub[i]].node[i];
AD[m][node_coarse] = AU[m][node_fine];
}
return;
}
void inject_vector(E,start_lev,AU,AD)
struct All_variables *E;
int start_lev;
double **AU,**AD; /* data on upper/lower mesh */
{
int i,j,m,el,node_coarse,node_fine,sl_minus,eqn_fine,eqn_coarse;
const int dims = E->mesh.nsd;
const int ends = enodes[dims];
void gather_to_1layer_id ();
if(start_lev == E->mesh.levmin) {
fprintf(E->fp,"Warning, attempting to project below lowest level\n");
return;
}
sl_minus = start_lev-1;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=E->lmesh.NEL[sl_minus];el++)
for(i=1;i<=ends;i++) {
node_coarse = E->IEN[sl_minus][m][el].node[i];
node_fine=E->IEN[start_lev][m][E->EL[sl_minus][m][el].sub[i]].node[i];
for (j=1;j<=dims;j++) {
eqn_fine = E->ID[start_lev][m][node_fine].doff[j];
eqn_coarse = E->ID[sl_minus][m][node_coarse].doff[j];
AD[m][eqn_coarse] = AU[m][eqn_fine];
}
}
return;
}
/* =====================================================
Function to inject data from coarse to fine grid (i.e.
just dropping values at shared grid points.
===================================================== */
void un_inject_vector(E,start_lev,AD,AU)
struct All_variables *E;
int start_lev;
double **AU,**AD; /* data on upper/lower mesh */
{
int i,m;
int el,node,node_plus;
int eqn1,eqn_plus1;
int eqn2,eqn_plus2;
int eqn3,eqn_plus3;
const int dims = E->mesh.nsd;
const int ends = enodes[dims];
const int sl_plus = start_lev+1;
const int neq = E->lmesh.NEQ[sl_plus];
const int nels = E->lmesh.NEL[start_lev];
assert(start_lev != E->mesh.levmax /* un_injection */);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<neq;i++)
AU[m][i]=0.0;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels;el++)
for(i=1;i<=ENODES3D;i++) {
node = E->IEN[start_lev][m][el].node[i];
node_plus=E->IEN[sl_plus][m][E->EL[start_lev][m][el].sub[i]].node[i];
eqn1 = E->ID[start_lev][m][node].doff[1];
eqn2 = E->ID[start_lev][m][node].doff[2];
eqn3 = E->ID[start_lev][m][node].doff[3];
eqn_plus1 = E->ID[sl_plus][m][node_plus].doff[1];
eqn_plus2 = E->ID[sl_plus][m][node_plus].doff[2];
eqn_plus3 = E->ID[sl_plus][m][node_plus].doff[3];
AU[m][eqn_plus1] = AD[m][eqn1];
AU[m][eqn_plus2] = AD[m][eqn2];
AU[m][eqn_plus3] = AD[m][eqn3];
}
return;
}
/* =======================================================================================
Interpolation from coarse grid to fine. See the appology attached to project() if you get
stressed out by node based assumptions. If it makes you feel any better, I don't like
it much either.
======================================================================================= */
void interp_vector(E,start_lev,AD,AU)
struct All_variables *E;
int start_lev;
double **AD,**AU; /* data on upper/lower mesh */
{
void un_inject_vector();
void fill_in_gaps();
void from_rtf_to_xyz();
void from_xyz_to_rtf();
void scatter_to_nlayer_id();
int i,j,k,m;
float x1,x2;
float n1,n2;
int noxz,node0,node1,node2;
int eqn0,eqn1,eqn2;
const int level = start_lev + 1;
const int dims =E->mesh.nsd;
const int ends= enodes[dims];
const int nox = E->lmesh.NOX[level];
const int noz = E->lmesh.NOZ[level];
const int noy = E->lmesh.NOY[level];
const int high_eqn = E->lmesh.NEQ[level];
if (start_lev==E->mesh.levmax) return;
from_rtf_to_xyz(E,start_lev,AD,AU); /* transform in xyz coordinates */
un_inject_vector(E,start_lev,AU,E->temp); /* information from lower level */
fill_in_gaps(E,E->temp,level);
from_xyz_to_rtf(E,level,E->temp,AU); /* get back to rtf coordinates */
return;
}
/* ==============================================
function to project viscosity down to all the
levels in the problem. (no gaps for vbcs)
============================================== */
void project_viscosity(E)
struct All_variables *E;
{
int lv,i,j,k,m,sl_minus;
void inject_scalar();
void project_scalar();
void project_scalar_e();
void inject_scalar_e();
void visc_from_gint_to_nodes();
void visc_from_nodes_to_gint();
void visc_from_gint_to_ele();
void visc_from_ele_to_gint();
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
float *viscU[NCS],*viscD[NCS];
lv = E->mesh.levmax;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
viscU[m]=(float *)malloc((1+E->lmesh.NNO[lv])*sizeof(float));
viscD[m]=(float *)malloc((1+vpts*E->lmesh.NEL[lv-1])*sizeof(float));
}
#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC /* allow for anisotropy */
if(E->viscosity.allow_anisotropic_viscosity){
for(lv=E->mesh.levmax;lv>E->mesh.levmin;lv--) {
sl_minus = lv -1;
if (E->viscosity.smooth_cycles==1) {
visc_from_gint_to_nodes(E,E->EVI[lv],viscU,lv); /* isotropic */
project_scalar(E,lv,viscU,viscD);
visc_from_nodes_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);
/* anisotropic */
visc_from_gint_to_nodes(E,E->EVI2[lv],viscU,lv);project_scalar(E,lv,viscU,viscD);visc_from_nodes_to_gint(E,viscD,E->EVI2[sl_minus],sl_minus);
visc_from_gint_to_nodes(E,E->EVIn1[lv],viscU,lv);project_scalar(E,lv,viscU,viscD);visc_from_nodes_to_gint(E,viscD,E->EVIn1[sl_minus],sl_minus);
visc_from_gint_to_nodes(E,E->EVIn2[lv],viscU,lv);project_scalar(E,lv,viscU,viscD);visc_from_nodes_to_gint(E,viscD,E->EVIn2[sl_minus],sl_minus);
visc_from_gint_to_nodes(E,E->EVIn3[lv],viscU,lv);project_scalar(E,lv,viscU,viscD);visc_from_nodes_to_gint(E,viscD,E->EVIn3[sl_minus],sl_minus);
}
else if (E->viscosity.smooth_cycles==2) {
visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv); /* isotropic */
inject_scalar_e(E,lv,viscU,E->EVI[sl_minus]);
/* anisotropic */
visc_from_gint_to_ele(E,E->EVI2[lv],viscU,lv);inject_scalar_e(E,lv,viscU,E->EVI2[sl_minus]);
visc_from_gint_to_ele(E,E->EVIn1[lv],viscU,lv);inject_scalar_e(E,lv,viscU,E->EVIn1[sl_minus]);
visc_from_gint_to_ele(E,E->EVIn2[lv],viscU,lv);inject_scalar_e(E,lv,viscU,E->EVIn2[sl_minus]);
visc_from_gint_to_ele(E,E->EVIn3[lv],viscU,lv);inject_scalar_e(E,lv,viscU,E->EVIn3[sl_minus]);
}
else if (E->viscosity.smooth_cycles==3) {
visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv);/* isotropic */
project_scalar_e(E,lv,viscU,viscD);
visc_from_ele_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);
/* anisotropic */
visc_from_gint_to_ele(E,E->EVI2[lv],viscU,lv);project_scalar_e(E,lv,viscU,viscD);visc_from_ele_to_gint(E,viscD,E->EVI2[sl_minus],sl_minus);
visc_from_gint_to_ele(E,E->EVIn1[lv],viscU,lv);project_scalar_e(E,lv,viscU,viscD);visc_from_ele_to_gint(E,viscD,E->EVIn1[sl_minus],sl_minus);
visc_from_gint_to_ele(E,E->EVIn2[lv],viscU,lv);project_scalar_e(E,lv,viscU,viscD);visc_from_ele_to_gint(E,viscD,E->EVIn2[sl_minus],sl_minus);
visc_from_gint_to_ele(E,E->EVIn3[lv],viscU,lv);project_scalar_e(E,lv,viscU,viscD);visc_from_ele_to_gint(E,viscD,E->EVIn3[sl_minus],sl_minus);
}
else if (E->viscosity.smooth_cycles==0) {
inject_scalar(E,lv,E->VI[lv],E->VI[sl_minus]);/* isotropic */
visc_from_nodes_to_gint(E,E->VI[sl_minus],E->EVI[sl_minus],sl_minus);
/* anisotropic */
inject_scalar(E,lv,E->VI2[lv],E->VI2[sl_minus]);visc_from_nodes_to_gint(E,E->VI2[sl_minus],E->EVI2[sl_minus],sl_minus);
inject_scalar(E,lv,E->VIn1[lv],E->VIn1[sl_minus]);visc_from_nodes_to_gint(E,E->VIn1[sl_minus],E->EVIn1[sl_minus],sl_minus);
inject_scalar(E,lv,E->VIn2[lv],E->VIn2[sl_minus]);visc_from_nodes_to_gint(E,E->VIn2[sl_minus],E->EVIn2[sl_minus],sl_minus);
inject_scalar(E,lv,E->VIn3[lv],E->VIn3[sl_minus]);visc_from_nodes_to_gint(E,E->VIn3[sl_minus],E->EVIn3[sl_minus],sl_minus);
}
normalize_director_at_gint(E,E->EVIn1[sl_minus],E->EVIn2[sl_minus],E->EVIn3[sl_minus],sl_minus);
}
}else{
#endif
/* regular operation, isotropic viscosity */
for(lv=E->mesh.levmax;lv>E->mesh.levmin;lv--) {
sl_minus = lv -1;
if (E->viscosity.smooth_cycles==1) {
visc_from_gint_to_nodes(E,E->EVI[lv],viscU,lv);
project_scalar(E,lv,viscU,viscD);
visc_from_nodes_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);
}
else if (E->viscosity.smooth_cycles==2) {
visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv);
inject_scalar_e(E,lv,viscU,E->EVI[sl_minus]);
}
else if (E->viscosity.smooth_cycles==3) {
visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv);
project_scalar_e(E,lv,viscU,viscD);
visc_from_ele_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);
}
else if (E->viscosity.smooth_cycles==0) {
/* visc_from_gint_to_nodes(E,E->EVI[lv],viscU,lv);
inject_scalar(E,lv,viscU,viscD);
visc_from_nodes_to_gint(E,viscD,E->EVI[sl_minus],sl_minus); */
inject_scalar(E,lv,E->VI[lv],E->VI[sl_minus]);
visc_from_nodes_to_gint(E,E->VI[sl_minus],E->EVI[sl_minus],sl_minus);
}
/* for(m=1;m<=E->sphere.caps_per_proc;m++) {
for (i=1;i<=E->lmesh.NEL[lv-1];i++)
fprintf (E->fp_out,"%d %g\n",i,viscD[m][i]);
for (i=1;i<=E->lmesh.NEL[lv];i++)
fprintf (E->fp_out,"%d %g\n",i,viscU[m][i]);
}
*/
}
#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
}
#endif
for(m=1;m<=E->sphere.caps_per_proc;m++) {
free((void *)viscU[m]);
free((void *)viscD[m]);
}
return;
}
/* ==================================================== */
void inject_scalar_e(E,start_lev,AU,AD)
struct All_variables *E;
int start_lev;
float **AU,**AD; /* data on upper/lower mesh */
{
int i,j,m;
int el,node,e;
float average,w;
void gather_to_1layer_ele ();
const int sl_minus = start_lev-1;
const int nels_minus=E->lmesh.NEL[start_lev-1];
const int dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
const int vpts=vpoints[E->mesh.nsd];
const int n_minus=nels_minus*vpts;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=n_minus;i++)
AD[m][i] = 0.0;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels_minus;el++)
for(i=1;i<=ENODES3D;i++) {
e = E->EL[sl_minus][m][el].sub[i];
AD[m][(el-1)*vpts+i] = AU[m][e];
}
return;
}
/* ==================================================== */
void project_scalar_e(E,start_lev,AU,AD)
struct All_variables *E;
int start_lev;
float **AU,**AD; /* data on upper/lower mesh */
{
int i,j,m;
int el,node,e;
float average,w;
void gather_to_1layer_ele ();
const int sl_minus = start_lev-1;
const int nels_minus=E->lmesh.NEL[start_lev-1];
const int dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
const double weight=(double) 1.0/ends;
const int vpts=vpoints[E->mesh.nsd];
const int n_minus=nels_minus*vpts;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=n_minus;i++)
AD[m][i] = 0.0;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels_minus;el++) {
average=0.0;
for(i=1;i<=ENODES3D;i++) {
e = E->EL[sl_minus][m][el].sub[i];
average += AU[m][e];
}
AD[m][el] = average*weight;
}
return;
}
/* ==================================================== */
void project_scalar(E,start_lev,AU,AD)
struct All_variables *E;
int start_lev;
float **AU,**AD; /* data on upper/lower mesh */
{
int i,j,m;
int el,node,node1;
float average,w;
void gather_to_1layer_node ();
const int sl_minus = start_lev-1;
const int nno_minus=E->lmesh.NNO[start_lev-1];
const int nels_minus=E->lmesh.NEL[start_lev-1];
const int dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
const double weight=(double) 1.0/ends;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nno_minus;i++)
AD[m][i] = 0.0;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels_minus;el++)
for(i=1;i<=ENODES3D;i++) {
average=0.0;
node1 = E->EL[sl_minus][m][el].sub[i];
for(j=1;j<=ENODES3D;j++) {
node=E->IEN[start_lev][m][node1].node[j];
average += AU[m][node];
}
w=weight*average;
node= E->IEN[sl_minus][m][el].node[i];
AD[m][node] += w * E->TWW[sl_minus][m][el].node[i];
}
(E->exchange_node_f)(E,AD,sl_minus);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nno_minus;i++) {
AD[m][i] *= E->MASS[sl_minus][m][i];
}
return;
}
/* this is prefered scheme with averages */
void project_vector(E,start_lev,AU,AD,ic)
struct All_variables *E;
int start_lev,ic;
double **AU,**AD; /* data on upper/lower mesh */
{
int i,j,m;
int el,node1,node,e1;
int eqn1,eqn_minus1;
int eqn2,eqn_minus2;
int eqn3,eqn_minus3;
double average1,average2,average3,w,weight;
float CPU_time(),time;
void from_rtf_to_xyz();
void from_xyz_to_rtf();
void gather_to_1layer_id ();
const int sl_minus = start_lev-1;
const int neq_minus=E->lmesh.NEQ[start_lev-1];
const int nno_minus=E->lmesh.NNO[start_lev-1];
const int nels_minus=E->lmesh.NEL[start_lev-1];
const int dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
if (ic==1)
weight = 1.0;
else
weight=(double) 1.0/ends;
if (start_lev==E->mesh.levmin) return;
/* convert into xyz coordinates */
from_rtf_to_xyz(E,start_lev,AU,E->temp);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=0;i<neq_minus;i++)
E->temp1[m][i] = 0.0;
/* smooth in xyz coordinates */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels_minus;el++)
for(i=1;i<=ENODES3D;i++) {
node= E->IEN[sl_minus][m][el].node[i];
average1=average2=average3=0.0;
e1 = E->EL[sl_minus][m][el].sub[i];
for(j=1;j<=ENODES3D;j++) {
node1=E->IEN[start_lev][m][e1].node[j];
average1 += E->temp[m][E->ID[start_lev][m][node1].doff[1]];
average2 += E->temp[m][E->ID[start_lev][m][node1].doff[2]];
average3 += E->temp[m][E->ID[start_lev][m][node1].doff[3]];
}
w = weight*E->TWW[sl_minus][m][el].node[i];
E->temp1[m][E->ID[sl_minus][m][node].doff[1]] += w * average1;
E->temp1[m][E->ID[sl_minus][m][node].doff[2]] += w * average2;
E->temp1[m][E->ID[sl_minus][m][node].doff[3]] += w * average3;
}
(E->solver.exchange_id_d)(E, E->temp1, sl_minus);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nno_minus;i++) {
E->temp1[m][E->ID[sl_minus][m][i].doff[1]] *= E->MASS[sl_minus][m][i];
E->temp1[m][E->ID[sl_minus][m][i].doff[2]] *= E->MASS[sl_minus][m][i];
E->temp1[m][E->ID[sl_minus][m][i].doff[3]] *= E->MASS[sl_minus][m][i];
}
/* back into rtf coordinates */
from_xyz_to_rtf(E,sl_minus,E->temp1,AD);
return;
}
/* ================================================= */
void from_xyz_to_rtf(E,level,xyz,rtf)
struct All_variables *E;
int level;
double **rtf,**xyz;
{
int i,j,m,eqn1,eqn2,eqn3;
double cost,cosf,sint,sinf;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for (i=1;i<=E->lmesh.NNO[level];i++) {
eqn1 = E->ID[level][m][i].doff[1];
eqn2 = E->ID[level][m][i].doff[2];
eqn3 = E->ID[level][m][i].doff[3];
sint = E->SinCos[level][m][0][i];
sinf = E->SinCos[level][m][1][i];
cost = E->SinCos[level][m][2][i];
cosf = E->SinCos[level][m][3][i];
rtf[m][eqn1] = xyz[m][eqn1]*cost*cosf
+ xyz[m][eqn2]*cost*sinf
- xyz[m][eqn3]*sint;
rtf[m][eqn2] = -xyz[m][eqn1]*sinf
+ xyz[m][eqn2]*cosf;
rtf[m][eqn3] = xyz[m][eqn1]*sint*cosf
+ xyz[m][eqn2]*sint*sinf
+ xyz[m][eqn3]*cost;
}
return;
}
/* ================================================= */
void from_rtf_to_xyz(E,level,rtf,xyz)
struct All_variables *E;
int level;
double **rtf,**xyz;
{
int i,j,m,eqn1,eqn2,eqn3;
double cost,cosf,sint,sinf;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for (i=1;i<=E->lmesh.NNO[level];i++) {
eqn1 = E->ID[level][m][i].doff[1];
eqn2 = E->ID[level][m][i].doff[2];
eqn3 = E->ID[level][m][i].doff[3];
sint = E->SinCos[level][m][0][i];
sinf = E->SinCos[level][m][1][i];
cost = E->SinCos[level][m][2][i];
cosf = E->SinCos[level][m][3][i];
xyz[m][eqn1] = rtf[m][eqn1]*cost*cosf
- rtf[m][eqn2]*sinf
+ rtf[m][eqn3]*sint*cosf;
xyz[m][eqn2] = rtf[m][eqn1]*cost*sinf
+ rtf[m][eqn2]*cosf
+ rtf[m][eqn3]*sint*sinf;
xyz[m][eqn3] = -rtf[m][eqn1]*sint
+ rtf[m][eqn3]*cost;
}
return;
}
/* ========================================================== */
void fill_in_gaps(E,temp,level)
struct All_variables *E;
int level;
double **temp;
{
int i,j,k,m;
float x1,x2;
float n1,n2;
int rnoz,noxz,node0,node1,node2;
int eqn0,eqn1,eqn2;
const int dims =E->mesh.nsd;
const int ends= enodes[dims];
const int nox = E->lmesh.NOX[level];
const int noz = E->lmesh.NOZ[level];
const int noy = E->lmesh.NOY[level];
const int sl_minus = level-1;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
n1 = n2 =0.5;
noxz = nox*noz;
for(k=1;k<=noy;k+=2) /* Fill in gaps in x direction */
for(j=1;j<=noz;j+=2)
for(i=2;i<nox;i+=2) {
node0 = j + (i-1)*noz + (k-1)*noxz; /* this node */
node1 = node0 - noz;
node2 = node0 + noz;
/* now for each direction */
eqn0=E->ID[level][m][node0].doff[1];
eqn1=E->ID[level][m][node1].doff[1];
eqn2=E->ID[level][m][node2].doff[1];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
eqn0=E->ID[level][m][node0].doff[2];
eqn1=E->ID[level][m][node1].doff[2];
eqn2=E->ID[level][m][node2].doff[2];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
eqn0=E->ID[level][m][node0].doff[3];
eqn1=E->ID[level][m][node1].doff[3];
eqn2=E->ID[level][m][node2].doff[3];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
}
n1 = n2 =0.5;
for(i=1;i<=nox;i++) /* Fill in gaps in y direction */
for(j=1;j<=noz;j+=2)
for(k=2;k<noy;k+=2) {
node0 = j + (i-1)*noz + (k-1)*noxz; /* this node */
node1 = node0 - noxz;
node2 = node0 + noxz;
eqn0=E->ID[level][m][node0].doff[1];
eqn1=E->ID[level][m][node1].doff[1];
eqn2=E->ID[level][m][node2].doff[1];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
eqn0=E->ID[level][m][node0].doff[2];
eqn1=E->ID[level][m][node1].doff[2];
eqn2=E->ID[level][m][node2].doff[2];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
eqn0=E->ID[level][m][node0].doff[3];
eqn1=E->ID[level][m][node1].doff[3];
eqn2=E->ID[level][m][node2].doff[3];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
}
for(j=2;j<noz;j+=2) {
x1 = E->sphere.R[level][j] - E->sphere.R[level][j-1];
x2 = E->sphere.R[level][j+1] - E->sphere.R[level][j];
n1 = x2/(x1+x2);
n2 = 1.0-n1;
for(k=1;k<=noy;k++) /* Fill in gaps in z direction */
for(i=1;i<=nox;i++) {
node0 = j + (i-1)*noz + (k-1)*noxz; /* this node */
node1 = node0 - 1;
node2 = node0 + 1;
eqn0=E->ID[level][m][node0].doff[1];
eqn1=E->ID[level][m][node1].doff[1];
eqn2=E->ID[level][m][node2].doff[1];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
eqn0=E->ID[level][m][node0].doff[2];
eqn1=E->ID[level][m][node1].doff[2];
eqn2=E->ID[level][m][node2].doff[2];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
eqn0=E->ID[level][m][node0].doff[3];
eqn1=E->ID[level][m][node1].doff[3];
eqn2=E->ID[level][m][node2].doff[3];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
}
}
} /* end for m */
return;
}
Computing file changes ...