[cig-commits] r18812 - mc/3D/CitcomS/branches/eheien_dev/lib
emheien at geodynamics.org
emheien at geodynamics.org
Tue Aug 2 09:52:33 PDT 2011
Author: emheien
Date: 2011-08-02 09:52:33 -0700 (Tue, 02 Aug 2011)
New Revision: 18812
Modified:
mc/3D/CitcomS/branches/eheien_dev/lib/Advection_diffusion.c
mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c
mc/3D/CitcomS/branches/eheien_dev/lib/Global_operations.c
mc/3D/CitcomS/branches/eheien_dev/lib/Nodal_mesh.c
mc/3D/CitcomS/branches/eheien_dev/lib/Process_buoyancy.c
mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c
mc/3D/CitcomS/branches/eheien_dev/lib/Topo_gravity.c
mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c
mc/3D/CitcomS/branches/eheien_dev/lib/Viscosity_structures.c
mc/3D/CitcomS/branches/eheien_dev/lib/global_defs.h
mc/3D/CitcomS/branches/eheien_dev/lib/prototypes.h
mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h
Log:
Further work in converting tracer code to C++
Added classes for triangular elements, UV coordinates
Simplified shape functions for tracers
Merged regional/full keep_within_bounds functions
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Advection_diffusion.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Advection_diffusion.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -143,8 +143,6 @@
{
int i,d,n,nel,el,node,m;
- void velo_from_element();
-
float adv_timestep;
float ts,uc1,uc2,uc3,uc,size,step,VV[4][9];
@@ -390,7 +388,6 @@
double diff, int bc, unsigned int **FLAGS)
{
void get_rtf_at_vpts();
- void velo_from_element();
int el,e,a,i,a1,m;
double Eres[9],rtf[4][9]; /* correction to the (scalar) Tdot field */
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -65,9 +65,6 @@
const CapBoundary bounds);
static double findradial(CartesianCoord vec,
BoundaryPoint bp);
-static void fix_radius(struct All_variables *E,
- SphericalCoord &sc,
- CartesianCoord &cc);
static void fix_angle(double *angle);
static int iget_radial_element(struct All_variables *E,
int j, int iel,
@@ -590,7 +587,7 @@
sc = new_tracer.getSphericalPos();
cc = new_tracer.getCartesianPos();
- iel=(E->trace.iget_element)(E,j,-99,cc,sc);
+ iel=(E->trace.iget_element)(E,j,UNDEFINED_ELEMENT,cc,sc);
if (iel<1) {
fprintf(E->trace.fpt,"Error(lost souls) - element not here?\n");
@@ -715,8 +712,8 @@
/* 5 6 5 7 */
/* 6 7 6 8 */
-void full_get_shape_functions(struct All_variables *E,
- double shp[9], ElementID nelem,
+void full_get_shape_functions(struct All_variables *E, int &shape_iwedge,
+ double shp[6], ElementID nelem,
SphericalCoord sc)
{
const int j = 1;
@@ -727,8 +724,8 @@
int itry;
CoordUV uv;
- double shape2d[4];
- double shaperad[3];
+ double shape2d[3];
+ double shaperad[2];
double shape[7];
int maxlevel=E->mesh.levmax;
@@ -746,7 +743,7 @@
/* Check first wedge (1 of 2) */
- iwedge=1;
+ iwedge=0;
next_wedge:
@@ -757,9 +754,9 @@
/* if any shape functions are negative, goto next wedge */
- if (shape2d[1]<eps||shape2d[2]<eps||shape2d[3]<eps)
+ if (shape2d[0]<eps||shape2d[1]<eps||shape2d[2]<eps)
{
- inum=inum+1;
+ inum++;
/* AKMA clean this up */
if (inum>3)
{
@@ -771,12 +768,12 @@
{
CartesianCoord cc;
fprintf(E->trace.fpt,"ERROR(gnomonic_interpolation)-inum>1\n");
- fprintf(E->trace.fpt,"shape %f %f %f\n",shape2d[1],shape2d[2],shape2d[3]);
+ fprintf(E->trace.fpt,"shape %f %f %f\n",shape2d[0],shape2d[1],shape2d[2]);
fprintf(E->trace.fpt,"u %f v %f element: %d \n",uv.u,uv.v, nelem);
fprintf(E->trace.fpt,"Element uv boundaries: \n");
for(kk=1;kk<=4;kk++) {
i = (E->ien[j][nelem].node[kk] - 1) / E->lmesh.noz + 1;
- fprintf(E->trace.fpt,"%d: U: %f V:%f\n",kk,E->gnomonic[i].u,E->gnomonic[i].v);
+ fprintf(E->trace.fpt,"%d: U: %f V:%f\n",kk,(*E->trace.gnomonic)[i].u,(*E->trace.gnomonic)[i].v);
}
fprintf(E->trace.fpt,"theta: %f phi: %f rad: %f\n",sc._theta,sc._phi,sc._rad);
fprintf(E->trace.fpt,"Element theta-phi boundaries: \n");
@@ -787,7 +784,7 @@
ival=icheck_element(E,j,nelem,cc,sc._rad);
fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
- ival=(E->trace.iget_element)(E, j, -99, cc, sc);
+ ival=(E->trace.iget_element)(E, j, UNDEFINED_ELEMENT, cc, sc);
fprintf(E->trace.fpt,"New Element?: %d\n",ival);
ival=icheck_column_neighbors(E,j,nelem,cc,sc._rad);
@@ -804,7 +801,7 @@
exit(10);
}
- iwedge=2;
+ iwedge=1;
goto next_wedge;
}
@@ -819,17 +816,17 @@
/* Sum of shape functions is 1 */
- shp[0]=iwedge;
- shp[1]=shaperad[1]*shape2d[1];
- shp[2]=shaperad[1]*shape2d[2];
- shp[3]=shaperad[1]*shape2d[3];
- shp[4]=shaperad[2]*shape2d[1];
- shp[5]=shaperad[2]*shape2d[2];
- shp[6]=shaperad[2]*shape2d[3];
+ shape_iwedge = iwedge;
+ shp[0]=shaperad[0]*shape2d[0];
+ shp[1]=shaperad[0]*shape2d[1];
+ shp[2]=shaperad[0]*shape2d[2];
+ shp[3]=shaperad[1]*shape2d[0];
+ shp[4]=shaperad[1]*shape2d[1];
+ shp[5]=shaperad[1]*shape2d[2];
/** debug **
fprintf(E->trace.fpt, "shp: %e %e %e %e %e %e\n",
- shp[1], shp[2], shp[3], shp[4], shp[5], shp[6]);
+ shp[0], shp[1], shp[2], shp[3], shp[4], shp[5]);
/**/
}
@@ -856,73 +853,40 @@
/* 6 7 6 8 */
CartesianCoord full_get_velocity(struct All_variables *E,
- int j, ElementID nelem,
- SphericalCoord sc)
+ int j, ElementID nelem,
+ SphericalCoord sc)
{
- int iwedge;
+ int iwedge, i;
const int sphere_key = 0;
-
- double shape[9];
- double VV[4][9];
- double vx[6],vy[6],vz[6];
-
- full_get_shape_functions(E, shape, nelem, sc);
- iwedge=shape[0];
-
+
+ double shape[6];
+ CartesianCoord VV[9], vel[6], result;
+
+ full_get_shape_functions(E, iwedge, shape, nelem, sc);
+
/* get cartesian velocity */
velo_from_element_d(E, VV, j, nelem, sphere_key);
-
+
/* depending on wedge, set up velocity points */
-
- if (iwedge==1)
- {
- vx[0]=VV[1][1];
- vx[1]=VV[1][2];
- vx[2]=VV[1][3];
- vx[3]=VV[1][5];
- vx[4]=VV[1][6];
- vx[5]=VV[1][7];
- vy[0]=VV[2][1];
- vy[1]=VV[2][2];
- vy[2]=VV[2][3];
- vy[3]=VV[2][5];
- vy[4]=VV[2][6];
- vy[5]=VV[2][7];
- vz[0]=VV[3][1];
- vz[1]=VV[3][2];
- vz[2]=VV[3][3];
- vz[3]=VV[3][5];
- vz[4]=VV[3][6];
- vz[5]=VV[3][7];
- }
- if (iwedge==2)
- {
- vx[0]=VV[1][1];
- vx[1]=VV[1][3];
- vx[2]=VV[1][4];
- vx[3]=VV[1][5];
- vx[4]=VV[1][7];
- vx[5]=VV[1][8];
- vy[0]=VV[2][1];
- vy[1]=VV[2][3];
- vy[2]=VV[2][4];
- vy[3]=VV[2][5];
- vy[4]=VV[2][7];
- vy[5]=VV[2][8];
- vz[0]=VV[3][1];
- vz[1]=VV[3][3];
- vz[2]=VV[3][4];
- vz[3]=VV[3][5];
- vz[4]=VV[3][7];
- vz[5]=VV[3][8];
- }
-
- return CartesianCoord(vx[0]*shape[1]+vx[1]*shape[2]+shape[3]*vx[2]+
- vx[3]*shape[4]+vx[4]*shape[5]+shape[6]*vx[5],
- vy[0]*shape[1]+vy[1]*shape[2]+shape[3]*vy[2]+
- vy[3]*shape[4]+vy[4]*shape[5]+shape[6]*vy[5],
- vz[0]*shape[1]+vz[1]*shape[2]+shape[3]*vz[2]+
- vz[3]*shape[4]+vz[4]*shape[5]+shape[6]*vz[5]);
+
+ if (iwedge==0) {
+ vel[0]=VV[1];
+ vel[1]=VV[2];
+ vel[2]=VV[3];
+ vel[3]=VV[5];
+ vel[4]=VV[6];
+ vel[5]=VV[7];
+ } else if (iwedge==1) {
+ vel[0]=VV[1];
+ vel[1]=VV[3];
+ vel[2]=VV[4];
+ vel[3]=VV[5];
+ vel[4]=VV[7];
+ vel[5]=VV[8];
+ }
+
+ for (i=0;i<6;++i) result += vel[i] * shape[i];
+ return result;
}
/***************************************************************/
@@ -938,38 +902,17 @@
CoordUV uv,
int iwedge, double * shape2d)
{
-
- double a0,a1,a2;
/* convert nelem to surface element number */
int n = (nelem - 1) / E->lmesh.elz + 1;
- /* shape function 1 */
+ /* shape functions */
+ shape2d[0] = E->trace.shape_coefs[j][iwedge][n]->applyShapeFunc(0, uv);
+ shape2d[1] = E->trace.shape_coefs[j][iwedge][n]->applyShapeFunc(1, uv);
+ shape2d[2] = E->trace.shape_coefs[j][iwedge][n]->applyShapeFunc(2, uv);
- a0=E->trace.shape_coefs[j][iwedge][1][n];
- a1=E->trace.shape_coefs[j][iwedge][2][n];
- a2=E->trace.shape_coefs[j][iwedge][3][n];
-
- shape2d[1]=a0+a1*uv.u+a2*uv.v;
-
- /* shape function 2 */
-
- a0=E->trace.shape_coefs[j][iwedge][4][n];
- a1=E->trace.shape_coefs[j][iwedge][5][n];
- a2=E->trace.shape_coefs[j][iwedge][6][n];
-
- shape2d[2]=a0+a1*uv.u+a2*uv.v;
-
- /* shape function 3 */
-
- a0=E->trace.shape_coefs[j][iwedge][7][n];
- a1=E->trace.shape_coefs[j][iwedge][8][n];
- a2=E->trace.shape_coefs[j][iwedge][9][n];
-
- shape2d[3]=a0+a1*uv.u+a2*uv.v;
-
/** debug **
fprintf(E->trace.fpt, "el=%d els=%d iwedge=%d shape=(%e %e %e)\n",
- nelem, n, iwedge, shape2d[1], shape2d[2], shape2d[3]);
+ nelem, n, iwedge, shape2d[0], shape2d[1], shape2d[2]);
/**/
}
@@ -985,7 +928,6 @@
int node1,node5;
double rad1,rad5,f1,f2,delrad;
-
const double eps=1e-6;
double top_bound=1.0+eps;
double bottom_bound=0.0-eps;
@@ -1004,30 +946,27 @@
/* Save a small amount of computation here */
/* because f1+f2=1, shapes can be switched */
/*
- shaperad[1]=1.0-f1=1.0-(1.0-f2)=f2;
- shaperad[2]=1.0-f2=1.0-(10-f1)=f1;
+ shaperad[0]=1.0-f1=1.0-(1.0-f2)=f2;
+ shaperad[1]=1.0-f2=1.0-(10-f1)=f1;
*/
- shaperad[1]=f2;
- shaperad[2]=f1;
+ shaperad[0]=f2;
+ shaperad[1]=f1;
/* Some error control */
- if (shaperad[1]>(top_bound)||shaperad[1]<(bottom_bound)||
- shaperad[2]>(top_bound)||shaperad[2]<(bottom_bound))
+ if (shaperad[0]>top_bound || shaperad[0]<bottom_bound ||
+ shaperad[1]>top_bound || shaperad[1]<bottom_bound)
{
fprintf(E->trace.fpt,"ERROR(get_radial_shape)\n");
+ fprintf(E->trace.fpt,"shaperad[0]: %f \n",shaperad[0]);
fprintf(E->trace.fpt,"shaperad[1]: %f \n",shaperad[1]);
- fprintf(E->trace.fpt,"shaperad[2]: %f \n",shaperad[2]);
fflush(E->trace.fpt);
exit(10);
}
}
-
-
-
/**************************************************************/
/* SPHERICAL TO UV */
/* */
@@ -1044,10 +983,10 @@
/* theta_f and phi_f are the reference points of the cap */
- phi_f = E->gnomonic_reference_phi;
+ phi_f = E->trace.gnomonic_reference_phi;
- cos_theta_f = E->gnomonic[0].u;
- sin_theta_f = E->gnomonic[0].v;
+ cos_theta_f = (*E->trace.gnomonic)[0].u;
+ sin_theta_f = (*E->trace.gnomonic)[0].v;
cost=cos(sc._theta);
/*
@@ -1240,7 +1179,7 @@
/* Allocate memory for regnodetoel */
/* Regtoel is an integer array which represents nodes on */
/* the regular mesh. Each node on the regular mesh contains */
- /* the real element value if one exists (-99 otherwise) */
+ /* the real element value if one exists (UNDEFINED_ELEMENT otherwise) */
@@ -1252,11 +1191,11 @@
}
- /* Initialize regnodetoel - reg elements not used =-99 */
+ /* Initialize regnodetoel - reg elements not used = UNDEFINED_ELEMENT */
for (kk=1;kk<=numregnodes;kk++)
{
- E->trace.regnodetoel[j][kk]=-99;
+ E->trace.regnodetoel[j][kk]=UNDEFINED_ELEMENT;
}
/* Begin Mapping (only need to use surface elements) */
@@ -1350,7 +1289,7 @@
for (kk=1;kk<=numregnodes;kk++)
{
- E->trace.regnodetoel[j][kk]=-99;
+ E->trace.regnodetoel[j][kk]=UNDEFINED_ELEMENT;
/* find theta and phi for a given regular node */
@@ -1452,7 +1391,7 @@
{
for (kk=1;kk<=numregnodes;kk++)
{
- if (E->trace.regnodetoel[j][kk]!=-99)
+ if (E->trace.regnodetoel[j][kk]!=UNDEFINED_ELEMENT)
{
if ( (E->trace.regnodetoel[j][kk]<1)||(E->trace.regnodetoel[j][kk]>E->lmesh.nel) )
{
@@ -1492,11 +1431,11 @@
/* a real element. The number of real elements a regular */
/* element touches (one of its nodes are in) is ichoice. */
/* Special ichoice notes: */
- /* ichoice=-1 all regular element nodes = -99 (no elements) */
+ /* ichoice=-1 all regular element nodes = UNDEFINED_ELEMENT (no elements) */
/* ichoice=0 all 4 corners within one element */
/* ichoice=1 one element choice (diff from ichoice 0 in */
/* that perhaps one reg node is in an element */
- /* and the rest are not (-99). */
+ /* and the rest are not (UNDEFINED_ELEMENT). */
/* ichoice>1 Multiple elements to check */
numregel= E->trace.numregel[j];
@@ -1686,7 +1625,6 @@
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
-
isum=0;
for (kk=0;kk<=4;kk++) isum=isum+istat_ichoice[j][kk];
fprintf(E->trace.fpt,"\n\nInformation regarding number of real elements per regular elements\n");
@@ -1775,7 +1713,7 @@
}
}
- return -99;
+ return UNDEFINED_ELEMENT;
}
@@ -1783,7 +1721,7 @@
/* */
/* This function check all columns until the proper one for */
/* a point (x,y,z) is found. The surface element is returned */
-/* else -99 if can't be found. */
+/* else UNDEFINED_ELEMENT if can't be found. */
static int icheck_all_columns(struct All_variables *E,
int j,
@@ -1801,7 +1739,7 @@
if (icheck==1) return nel;
}
- return -99;
+ return UNDEFINED_ELEMENT;
}
@@ -2040,44 +1978,6 @@
}
-/************ FIX RADIUS ********************************************/
-/* This function moves particles back in bounds if they left */
-/* during advection */
-
-static void fix_radius(struct All_variables *E, SphericalCoord &sc, CartesianCoord &cc)
-{
- double sint,cost,sinf,cosf,rad;
- double max_radius, min_radius;
- int changed = 0;
-
- max_radius = E->sphere.ro - E->trace.box_cushion;
- min_radius = E->sphere.ri + E->trace.box_cushion;
-
- if (sc._rad > max_radius) {
- changed = 1;
- sc._rad=max_radius;
- rad=max_radius;
- }
- if (sc._rad < min_radius) {
- changed = 1;
- sc._rad=min_radius;
- rad=min_radius;
- }
-
- if (changed) {
- cost=cos(sc._theta);
- sint=sqrt(1.0-cost*cost);
- cosf=cos(sc._phi);
- sinf=sin(sc._phi);
- cc._x=rad*sint*cosf;
- cc._y=rad*sint*sinf;
- cc._z=rad*cost;
- }
-
- return;
-}
-
-
/******************************************************************/
/* FIX ANGLE */
/* */
@@ -2097,7 +1997,7 @@
/********** IGET ELEMENT *****************************************/
/* */
/* This function returns the the real element for a given point. */
-/* Returns -99 if not in this cap. */
+/* Returns UNDEFINED_ELEMENT if not in this cap. */
/* Returns -1 if in this cap but cannot find the element. */
/* iprevious_element, if known, is the last known element. If */
/* it is not known, input a negative number. */
@@ -2130,7 +2030,7 @@
if (E->parallel.nprocz>1)
{
ival=icheck_processor_shell(E,j,sc._rad);
- if (ival!=1) return -99;
+ if (ival!=1) return UNDEFINED_ELEMENT;
}
/* do quick search to see if element can be easily found. */
@@ -2143,7 +2043,7 @@
iregel=iget_regel(E,j,sc._theta,sc._phi,&ntheta,&nphi);
if (iregel<=0)
{
- return -99;
+ return UNDEFINED_ELEMENT;
}
@@ -2207,7 +2107,7 @@
/* check if still in cap */
ival=full_icheck_cap(E,0,cc,sc._rad);
- if (ival==0) return -99;
+ if (ival==0) return UNDEFINED_ELEMENT;
/* if here, still in cap (hopefully, without a doubt) */
@@ -2347,7 +2247,7 @@
/*********** IGET REGEL ******************************************/
/* */
/* This function returns the regular element in which a point */
-/* exists. If not found, returns -99. */
+/* exists. If not found, returns UNDEFINED_ELEMENT. */
/* npi and ntheta are modified for later use */
static int iget_regel(struct All_variables *E, int j,
@@ -2360,8 +2260,8 @@
/* first check whether theta is in range */
- if (theta<E->trace.thetamin[j]) return -99;
- if (theta>E->trace.thetamax[j]) return -99;
+ if (theta<E->trace.thetamin[j]) return UNDEFINED_ELEMENT;
+ if (theta>E->trace.thetamax[j]) return UNDEFINED_ELEMENT;
/* get ntheta, nphi on regular mesh */
@@ -2377,8 +2277,8 @@
/* check range to be sure */
- if (iregel>E->trace.numregel[j]) return -99;
- if (iregel<1) return -99;
+ if (iregel>E->trace.numregel[j]) return UNDEFINED_ELEMENT;
+ if (iregel<1) return UNDEFINED_ELEMENT;
return iregel;
}
@@ -2404,12 +2304,6 @@
double u, v, cosc, theta_f, phi_f, dphi, cosd;
double *cost, *sint, *cosf, *sinf;
- if ((E->gnomonic = (CITCOM_GNOMONIC*)malloc((E->lmesh.nsf+1)*sizeof(struct CITCOM_GNOMONIC)))
- == NULL) {
- fprintf(stderr,"Error(define uv)-not enough memory(a)\n");
- exit(10);
- }
-
sint = E->SinCos[lev][j][0];
sinf = E->SinCos[lev][j][1];
cost = E->SinCos[lev][j][2];
@@ -2419,7 +2313,7 @@
/* use the point at middle of the cap */
refnode = 1 + E->lmesh.noz * ((E->lmesh.noy / 2) * E->lmesh.nox
+ E->lmesh.nox / 2);
- phi_f = E->gnomonic_reference_phi = E->sx[j][2][refnode];
+ phi_f = E->trace.gnomonic_reference_phi = E->sx[j][2][refnode];
/** debug **
theta_f = E->sx[j][1][refnode];
@@ -2432,10 +2326,8 @@
/**/
/* store cos(theta_f) and sin(theta_f) */
- E->gnomonic[0].u = cost[refnode];
- E->gnomonic[0].v = sint[refnode];
+ E->trace.gnomonic->insert(std::make_pair(0, CoordUV(cost[refnode], sint[refnode])));
-
/* convert each nodal point to u and v */
for (i=1, n=1; i<=E->lmesh.nsf; i++, n+=E->lmesh.noz) {
@@ -2446,19 +2338,15 @@
v = (sint[refnode] * cost[n] - cost[refnode] * sint[n] * cosd)
/ cosc;
- E->gnomonic[i].u = u;
- E->gnomonic[i].v = v;
+ E->trace.gnomonic->insert(std::make_pair(i, CoordUV(u, v)));
/** debug **
fprintf(E->trace.fpt, "n=%d ns=%d cosc=%e (%e %e) -> (%e %e)\n",
n, i, cosc, E->sx[j][1][n], E->sx[j][2][n], u, v);
/**/
}
-
- return;
}
-
/***************************************************************/
/* DETERMINE SHAPE COEFFICIENTS */
/* */
@@ -2476,26 +2364,17 @@
int nelem, iwedge, kk, i;
int snode;
- double u[5], v[5];
- double x1 = 0.0;
- double x2 = 0.0;
- double x3 = 0.0;
- double y1 = 0.0;
- double y2 = 0.0;
- double y3 = 0.0;
- double delta, a0, a1, a2;
+ CoordUV uv[5], xy1, xy2, xy3;
/* first, allocate memory */
- for(iwedge=1; iwedge<=2; iwedge++) {
- for (kk=1; kk<=9; kk++) {
- if ((E->trace.shape_coefs[j][iwedge][kk] =
- (double *)malloc((E->lmesh.snel+1)*sizeof(double))) == NULL) {
- fprintf(E->trace.fpt,"ERROR(find shape coefs)-not enough memory(a)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ for(iwedge=0; iwedge<2; iwedge++) {
+ if ((E->trace.shape_coefs[j][iwedge] =
+ (TriElemLinearShapeFunc **)malloc((E->lmesh.snel+1)*sizeof(TriElemLinearShapeFunc*))) == NULL) {
+ fprintf(E->trace.fpt,"ERROR(find shape coefs)-not enough memory(a)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
}
for (i=1, nelem=1; i<=E->lmesh.snel; i++, nelem+=E->lmesh.elz) {
@@ -2504,65 +2383,27 @@
for(kk=1; kk<=4; kk++) {
snode = (E->ien[j][nelem].node[kk]-1) / E->lmesh.noz + 1;
- u[kk] = E->gnomonic[snode].u;
- v[kk] = E->gnomonic[snode].v;
+ uv[kk] = (*E->trace.gnomonic)[snode];
}
- for(iwedge=1; iwedge<=2; iwedge++) {
+ for(iwedge=0; iwedge<2; iwedge++) {
- if (iwedge == 1) {
- x1 = u[1];
- x2 = u[2];
- x3 = u[3];
- y1 = v[1];
- y2 = v[2];
- y3 = v[3];
+ if (iwedge == 0) {
+ xy1 = uv[1];
+ xy2 = uv[2];
+ xy3 = uv[3];
+ } else if (iwedge == 1) {
+ xy1 = uv[1];
+ xy2 = uv[3];
+ xy3 = uv[4];
}
- if (iwedge == 2) {
- x1 = u[1];
- x2 = u[3];
- x3 = u[4];
- y1 = v[1];
- y2 = v[3];
- y3 = v[4];
- }
+
+ E->trace.shape_coefs[j][iwedge][i] = new TriElemLinearShapeFunc(xy1, xy2, xy3);
- /* shape function 1 */
-
- delta = (x3-x2)*(y1-y2)-(y2-y3)*(x2-x1);
- a0 = (x2*y3-x3*y2)/delta;
- a1 = (y2-y3)/delta;
- a2 = (x3-x2)/delta;
-
- E->trace.shape_coefs[j][iwedge][1][i] = a0;
- E->trace.shape_coefs[j][iwedge][2][i] = a1;
- E->trace.shape_coefs[j][iwedge][3][i] = a2;
-
- /* shape function 2 */
-
- delta = (x3-x1)*(y2-y1)-(y1-y3)*(x1-x2);
- a0 = (x1*y3-x3*y1)/delta;
- a1 = (y1-y3)/delta;
- a2 = (x3-x1)/delta;
-
- E->trace.shape_coefs[j][iwedge][4][i] = a0;
- E->trace.shape_coefs[j][iwedge][5][i] = a1;
- E->trace.shape_coefs[j][iwedge][6][i] = a2;
-
- /* shape function 3 */
-
- delta = (x1-x2)*(y3-y2)-(y2-y1)*(x2-x3);
- a0 = (x2*y1-x1*y2)/delta;
- a1 = (y2-y1)/delta;
- a2 = (x1-x2)/delta;
-
- E->trace.shape_coefs[j][iwedge][7][i] = a0;
- E->trace.shape_coefs[j][iwedge][8][i] = a1;
- E->trace.shape_coefs[j][iwedge][9][i] = a2;
-
/** debug **
fprintf(E->trace.fpt, "el=%d els=%d iwedge=%d shape=(%e %e %e, %e %e %e, %e %e %e)\n",
nelem, i, iwedge,
+ E->trace.shape_coefs[j][iwedge][0][i],
E->trace.shape_coefs[j][iwedge][1][i],
E->trace.shape_coefs[j][iwedge][2][i],
E->trace.shape_coefs[j][iwedge][3][i],
@@ -2570,8 +2411,7 @@
E->trace.shape_coefs[j][iwedge][5][i],
E->trace.shape_coefs[j][iwedge][6][i],
E->trace.shape_coefs[j][iwedge][7][i],
- E->trace.shape_coefs[j][iwedge][8][i],
- E->trace.shape_coefs[j][iwedge][9][i]);
+ E->trace.shape_coefs[j][iwedge][8][i]);
/**/
} /* end wedge */
@@ -2579,20 +2419,6 @@
}
-/*********** KEEP WITHIN BOUNDS *****************************************/
-/* */
-/* This function makes sure the particle is within the sphere, and */
-/* phi and theta are within the proper degree range. */
-
-void full_keep_within_bounds(struct All_variables *E,
- CartesianCoord &cc,
- SphericalCoord &sc)
-{
- sc.constrainThetaPhi();
- fix_radius(E,sc,cc);
-}
-
-
/* &&&&&&&&&&&&&&&&&&&& ANALYTICAL TESTS &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&**************/
/**************** ANALYTICAL TEST *********************************************************/
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Global_operations.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Global_operations.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -845,7 +845,7 @@
const int ppts = PPOINTS3D;
const int vpts = VPOINTS3D;
const int sphere_key = 1;
- double VV[4][9];
+ CartesianCoord VV[9];
double rot, fr, tr;
double tmp, moment_of_inertia, rho;
@@ -891,8 +891,8 @@
}
for (j=1;j<=ppts;j++) {
for (i=1;i<=ends;i++) {
- vx[j] += VV[1][i]*E->N.ppt[GNPINDEX(i,j)];
- vy[j] += VV[2][i]*E->N.ppt[GNPINDEX(i,j)];
+ vx[j] += VV[i]._x*E->N.ppt[GNPINDEX(i,j)];
+ vy[j] += VV[i]._y*E->N.ppt[GNPINDEX(i,j)];
}
}
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Nodal_mesh.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Nodal_mesh.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Nodal_mesh.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -153,7 +153,7 @@
}
/* double prec version */
-void velo_from_element_d(struct All_variables *E, double VV[4][9], int m, int el, int sphere_key)
+void velo_from_element_d(struct All_variables *E, CartesianCoord VV[9], int m, int el, int sphere_key)
{
int a, node;
@@ -166,9 +166,9 @@
if (sphere_key)
for(a=1;a<=ends;a++) {
node = E->ien[m][el].node[a];
- VV[1][a] = E->sphere.cap[m].V[1][node];
- VV[2][a] = E->sphere.cap[m].V[2][node];
- VV[3][a] = E->sphere.cap[m].V[3][node];
+ VV[a] = CartesianCoord(E->sphere.cap[m].V[1][node],
+ E->sphere.cap[m].V[2][node],
+ E->sphere.cap[m].V[3][node]);
}
else {
for(a=1;a<=ends;a++) {
@@ -179,17 +179,16 @@
cost = E->SinCos[lev][m][2][node];
cosf = E->SinCos[lev][m][3][node];
- VV[1][a] = E->sphere.cap[m].V[1][node]*cost*cosf
- - E->sphere.cap[m].V[2][node]*sinf
- + E->sphere.cap[m].V[3][node]*sint*cosf;
- VV[2][a] = E->sphere.cap[m].V[1][node]*cost*sinf
- + E->sphere.cap[m].V[2][node]*cosf
- + E->sphere.cap[m].V[3][node]*sint*sinf;
- VV[3][a] = -E->sphere.cap[m].V[1][node]*sint
- + E->sphere.cap[m].V[3][node]*cost;
+ VV[a] = CartesianCoord(E->sphere.cap[m].V[1][node]*cost*cosf
+ - E->sphere.cap[m].V[2][node]*sinf
+ + E->sphere.cap[m].V[3][node]*sint*cosf,
+ E->sphere.cap[m].V[1][node]*cost*sinf
+ + E->sphere.cap[m].V[2][node]*cosf
+ + E->sphere.cap[m].V[3][node]*sint*sinf,
+ -E->sphere.cap[m].V[1][node]*sint
+ + E->sphere.cap[m].V[3][node]*cost);
}
}
- return;
}
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Process_buoyancy.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Process_buoyancy.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Process_buoyancy.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -57,7 +57,6 @@
float VV[4][9],u[9],T[9],dTdz[9],rho[9],area,uT;
float *sum_h;
- void velo_from_element();
void sum_across_surface();
void return_horiz_ave();
void return_horiz_ave_f();
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -126,7 +126,7 @@
/********** IGET ELEMENT *****************************************/
/* */
/* This function returns the the real element for a given point. */
-/* Returns -99 in not in this cap. */
+/* Returns UNDEFINED_ELEMENT in not in this cap. */
/* iprevious_element, if known, is the last known element. If */
/* it is not known, input a negative number. */
@@ -167,7 +167,7 @@
kk = isearch_all(E->trace.z_space, elz+1, sc._rad);
if (ii<0 || jj<0 || kk<0)
- return -99;
+ return UNDEFINED_ELEMENT;
return jj*elx*elz + ii*elz + kk + 1;
}
@@ -331,10 +331,11 @@
int m, int nelem,
SphericalCoord sc)
{
- double shp[8], VV[4][9], tmp;
- int n, d, node;
+ double shp[8], tmp;
+ CartesianCoord VV[9];
+ int n;
const int sphere_key = 0;
- double velocity_vector[3];
+ CartesianCoord velocity_vector;
/* get shape functions at (theta, phi, rad) */
regional_get_shape_functions(E, shp, nelem, sc);
@@ -348,17 +349,13 @@
Interpolate the velocity on the tracer position
***/
- for(d=0; d<3; d++)
- velocity_vector[d] = 0;
+ for(n=1; n<9; n++) {
+ velocity_vector += VV[n] * shp[n];
+ }
- for(d=0; d<3; d++) {
- for(n=0; n<8; n++)
- velocity_vector[d] += VV[d+1][n+1] * shp[n];
- }
-
/** debug **
- for(d=1; d<=3; d++) {
+ for(d=0; d<3; d++) {
fprintf(E->trace.fpt, "VV: %e %e %e %e %e %e %e %e: %e\n",
VV[d][1], VV[d][2], VV[d][3], VV[d][4],
VV[d][5], VV[d][6], VV[d][7], VV[d][8],
@@ -374,51 +371,10 @@
fflush(E->trace.fpt);
/**/
- return CartesianCoord(velocity_vector[0], velocity_vector[1], velocity_vector[2]);
+ return velocity_vector;
}
-void regional_keep_within_bounds(struct All_variables *E,
- CartesianCoord &cc,
- SphericalCoord &sc)
-{
- int changed = 0;
-
- if (sc._theta > E->control.theta_max - E->trace.box_cushion) {
- sc._theta = E->control.theta_max - E->trace.box_cushion;
- changed = 1;
- }
-
- if (sc._theta < E->control.theta_min + E->trace.box_cushion) {
- sc._theta = E->control.theta_min + E->trace.box_cushion;
- changed = 1;
- }
-
- if (sc._phi > E->control.fi_max - E->trace.box_cushion) {
- sc._phi = E->control.fi_max - E->trace.box_cushion;
- changed = 1;
- }
-
- if (sc._phi < E->control.fi_min + E->trace.box_cushion) {
- sc._phi = E->control.fi_min + E->trace.box_cushion;
- changed = 1;
- }
-
- if (sc._rad > E->sphere.ro - E->trace.box_cushion) {
- sc._rad = E->sphere.ro - E->trace.box_cushion;
- changed = 1;
- }
-
- if (sc._rad < E->sphere.ri + E->trace.box_cushion) {
- sc._rad = E->sphere.ri + E->trace.box_cushion;
- changed = 1;
- }
-
- if (changed)
- cc = sc.toCartesian();
-}
-
-
void regional_lost_souls(struct All_variables *E)
{
/* This part only works if E->sphere.caps_per_proc==1 */
@@ -762,7 +718,7 @@
new_tracer.readFromMem(&recv[ipos]);
/* found the element */
- iel = regional_iget_element(E, j, -99, CartesianCoord(0, 0, 0), sc);
+ iel = regional_iget_element(E, j, UNDEFINED_ELEMENT, CartesianCoord(0, 0, 0), sc);
if (iel<1) {
fprintf(E->trace.fpt, "Error(regional lost souls) - "
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Topo_gravity.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Topo_gravity.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -212,7 +212,6 @@
float** divv, float** vorv)
{
void get_rtf_at_vpts();
- void velo_from_element();
void stress_conform_bcs();
void construct_c3x3matrix_el();
void get_ba();
@@ -991,7 +990,6 @@
void get_global_1d_shape_fn_L();
void full_exchange_snode_f();
void regional_exchange_snode_f();
- void velo_from_element();
int a,address,el,elb,els,node,nodeb,nodes,i,j,k,l,m,n,count;
int nodel,nodem,nodesl,nodesm,nnsf,nel2;
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -183,32 +183,64 @@
return CartesianCoord(x,y,z);
}
+// Add a coordinate to this one by summing the components
+CartesianCoord & CartesianCoord::operator+=(const CartesianCoord &other) {
+ this->_x += other._x;
+ this->_y += other._y;
+ this->_z += other._z;
+ return *this;
+}
+
+// Get the difference of a coordinate from this one as vectors by subtracting components
+CartesianCoord & CartesianCoord::operator-=(const CartesianCoord &other) {
+ this->_x -= other._x;
+ this->_y -= other._y;
+ this->_z -= other._z;
+ return *this;
+}
+
+// Multiply each component by a constant factor
+CartesianCoord & CartesianCoord::operator*=(const double &val) {
+ this->_x *= val;
+ this->_y *= val;
+ this->_z *= val;
+ return *this;
+}
+
+// Divide each element by a constant factor
+CartesianCoord & CartesianCoord::operator/=(const double &val) {
+ this->_x /= val;
+ this->_y /= val;
+ this->_z /= val;
+ return *this;
+}
+
// Add two coordinates together as vectors by summing each of their components
const CartesianCoord CartesianCoord::operator+(const CartesianCoord &other) const {
- return CartesianCoord( this->_x+other._x,
- this->_y+other._y,
- this->_z+other._z);
+ CartesianCoord new_coord = *this;
+ new_coord += other;
+ return new_coord;
}
// Get the difference of two coordinates as vectors by subtracting each of their components
const CartesianCoord CartesianCoord::operator-(const CartesianCoord &other) const {
- return CartesianCoord( this->_x-other._x,
- this->_y-other._y,
- this->_z-other._z);
+ CartesianCoord new_coord = *this;
+ new_coord -= other;
+ return new_coord;
}
// Multiply each component by a constant factor
const CartesianCoord CartesianCoord::operator*(const double &val) const {
- return CartesianCoord( this->_x*val,
- this->_y*val,
- this->_z*val);
+ CartesianCoord new_coord = *this;
+ new_coord *= val;
+ return new_coord;
}
// Divide each element by a constant factor
const CartesianCoord CartesianCoord::operator/(const double &val) const {
- return CartesianCoord( this->_x/val,
- this->_y/val,
- this->_z/val);
+ CartesianCoord new_coord = *this;
+ new_coord /= val;
+ return new_coord;
}
// Return the cross product of this vector with another vector (b)
@@ -258,6 +290,46 @@
}
+TriElemLinearShapeFunc::TriElemLinearShapeFunc(CoordUV xy1, CoordUV xy2, CoordUV xy3) {
+ double delta, a0, a1, a2;
+
+ /* shape function 1 */
+ delta = (xy3.u-xy2.u)*(xy1.v-xy2.v)-(xy2.v-xy3.v)*(xy2.u-xy1.u);
+ a0 = (xy2.u*xy3.v-xy3.u*xy2.v)/delta;
+ a1 = (xy2.v-xy3.v)/delta;
+ a2 = (xy3.u-xy2.u)/delta;
+
+ sf[0][0] = a0;
+ sf[0][1] = a1;
+ sf[0][2] = a2;
+
+ /* shape function 2 */
+ delta = (xy3.u-xy1.u)*(xy2.v-xy1.v)-(xy1.v-xy3.v)*(xy1.u-xy2.u);
+ a0 = (xy1.u*xy3.v-xy3.u*xy1.v)/delta;
+ a1 = (xy1.v-xy3.v)/delta;
+ a2 = (xy3.u-xy1.u)/delta;
+
+ sf[1][0] = a0;
+ sf[1][1] = a1;
+ sf[1][2] = a2;
+
+ /* shape function 3 */
+
+ delta = (xy1.u-xy2.u)*(xy3.v-xy2.v)-(xy2.v-xy1.v)*(xy2.u-xy3.u);
+ a0 = (xy2.u*xy1.v-xy1.u*xy2.v)/delta;
+ a1 = (xy2.v-xy1.v)/delta;
+ a2 = (xy1.u-xy2.u)/delta;
+
+ sf[2][0] = a0;
+ sf[2][1] = a1;
+ sf[2][2] = a2;
+}
+
+double TriElemLinearShapeFunc::applyShapeFunc(const unsigned int i, const CoordUV uv) const {
+ assert(i < 3);
+ return sf[i][0]+sf[i][1]*uv.u+sf[i][2]*uv.v;
+}
+
void tracer_input(struct All_variables *E)
{
char message[100];
@@ -378,22 +450,21 @@
// a regional or full mantle simulation
void tracer_initial_settings(struct All_variables *E)
{
+ // Initialize CPU time counters
E->trace.advection_time = 0;
E->trace.find_tracers_time = 0;
E->trace.lost_souls_time = 0;
- if(E->parallel.nprocxy == 1) {
- E->trace.full_tracers = false;
-
- E->trace.keep_within_bounds = regional_keep_within_bounds;
+ // If we have more than 1 processor, we are in full mantle tracer mode
+ E->trace.full_tracers = (E->parallel.nprocxy != 1);
+
+ // Set up function pointers depending on which mode we are in
+ if(E->trace.full_tracers) {
+ E->trace.get_velocity = full_get_velocity;
+ E->trace.iget_element = full_iget_element;
+ } else {
E->trace.get_velocity = regional_get_velocity;
E->trace.iget_element = regional_iget_element;
- } else {
- E->trace.full_tracers = true;
-
- E->trace.keep_within_bounds = full_keep_within_bounds;
- E->trace.get_velocity = full_get_velocity;
- E->trace.iget_element = full_iget_element;
}
}
@@ -431,6 +502,9 @@
write_trace_instructions(E);
if (E->trace.full_tracers) {
+ // Initialize gnomonic node->coordinate mapping
+ E->trace.gnomonic = new std::map<int, CoordUV>;
+
/* Gnometric projection for velocity interpolation */
define_uv_space(E);
determine_shape_coefficients(E);
@@ -653,6 +727,83 @@
}
+/*********** KEEP WITHIN BOUNDS *****************************************/
+/* */
+/* This function makes sure the particle is within the sphere, and */
+/* phi and theta are within the proper degree range. */
+
+void keep_within_bounds(struct All_variables *E,
+ CartesianCoord &cc,
+ SphericalCoord &sc)
+{
+ double max_theta, min_theta;
+ double max_fi, min_fi;
+ double max_radius, min_radius;
+ double sint,cost,sinf,cosf;
+ int changed = 0;
+
+ // TODO: simplify and unify (if possible) the logic here
+ // If we're in full tracer mode, constrain phi and theta
+ if (E->trace.full_tracers) {
+ sc.constrainThetaPhi();
+ } else { // otherwise use the min/max values
+ max_theta = E->control.theta_max - E->trace.box_cushion;
+ min_theta = E->control.theta_min + E->trace.box_cushion;
+
+ max_fi = E->control.fi_max - E->trace.box_cushion;
+ min_fi = E->control.fi_min + E->trace.box_cushion;
+
+ if (sc._theta > max_theta) {
+ sc._theta = max_theta;
+ changed = 1;
+ }
+
+ if (sc._theta < min_theta) {
+ sc._theta = min_theta;
+ changed = 1;
+ }
+
+ if (sc._phi > max_fi) {
+ sc._phi = max_fi;
+ changed = 1;
+ }
+
+ if (sc._phi < min_fi) {
+ sc._phi = min_fi;
+ changed = 1;
+ }
+ }
+
+ // In all modes we are constrained by radius
+ max_radius = E->sphere.ro - E->trace.box_cushion;
+ min_radius = E->sphere.ri + E->trace.box_cushion;
+
+ if (sc._rad > max_radius) {
+ changed = 1;
+ sc._rad=max_radius;
+ }
+ if (sc._rad < min_radius) {
+ changed = 1;
+ sc._rad=min_radius;
+ }
+
+ // If the particle was moved, recalculate the Cartesian position
+ if (changed) {
+ if (E->trace.full_tracers) {
+ cost=cos(sc._theta);
+ sint=sqrt(1.0-cost*cost);
+ cosf=cos(sc._phi);
+ sinf=sin(sc._phi);
+ cc._x=sc._rad*sint*cosf;
+ cc._y=sc._rad*sint*sinf;
+ cc._z=sc._rad*cost;
+ } else {
+ cc = sc.toCartesian();
+ }
+ }
+}
+
+
/*********** PREDICT TRACERS **********************************************/
/* */
/* This function predicts tracers performing an euler step */
@@ -686,7 +837,7 @@
// Ensure tracer remains in the boundaries
sc_pred = cc_pred.toSpherical();
- (E->trace.keep_within_bounds)(E, cc_pred, sc_pred);
+ keep_within_bounds(E, cc_pred, sc_pred);
// Set new tracer coordinates
tr->setCoords(sc_pred, cc_pred);
@@ -735,7 +886,7 @@
new_coord_sph = new_coord.toSpherical();
// Ensure tracer remains in boundaries
- (E->trace.keep_within_bounds)(E, new_coord, new_coord_sph);
+ keep_within_bounds(E, new_coord, new_coord_sph);
// Fill in current positions (original values are no longer important)
tr->setCoords(new_coord_sph, new_coord);
@@ -752,7 +903,7 @@
/* */
/* This function finds tracer elements and moves tracers to */
/* other processor domains if necessary. */
-/* Array ielement is filled with elemental values. */
+/* Array ielement is filled with elemental values. */
static void find_tracers(struct All_variables *E)
{
@@ -790,7 +941,7 @@
tr->set_ielement(iel);
- if (iel == -99) {
+ if (iel == UNDEFINED_ELEMENT) {
/* tracer is inside other processors */
E->trace.escaped_tracers[j].push_back(*tr);
tr=E->trace.tracers[j].erase(tr);
@@ -866,20 +1017,27 @@
void initialize_tracers(struct All_variables *E)
{
+ // Initialize per-cap arrays for tracers and escaped tracers
E->trace.tracers = new TracerList[13];
E->trace.escaped_tracers = new TracerList[13];
- if (E->trace.ic_method==0)
- make_tracer_array(E);
- else if (E->trace.ic_method==1)
- read_tracer_file(E);
- else if (E->trace.ic_method==2)
- read_old_tracer_file(E);
- else {
- fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
+ // Call the appropriate function depending on the tracer input method
+ switch (E->trace.ic_method) {
+ case 0:
+ make_tracer_array(E);
+ break;
+ case 1:
+ read_tracer_file(E);
+ break;
+ case 2:
+ read_old_tracer_file(E);
+ break;
+ default:
+ fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ break;
+ }
/* total number of tracers */
@@ -998,16 +1156,16 @@
if (new_sc._rad<E->sx[j][3][1]) continue;
/* check if in current cap */
- if (E->parallel.nprocxy==1)
+ if (E->trace.full_tracers)
+ ival=full_icheck_cap(E,0,new_cc,new_sc._rad);
+ else
ival=regional_icheck_cap(E,0,new_sc,new_sc._rad);
- else
- ival=full_icheck_cap(E,0,new_cc,new_sc._rad);
if (ival!=1) continue;
/* Made it, so record tracer information */
- (E->trace.keep_within_bounds)(E, new_cc, new_sc);
+ keep_within_bounds(E, new_cc, new_sc);
Tracer new_tracer;
@@ -1089,7 +1247,7 @@
/* make sure theta, phi is in range, and radius is within bounds */
- (E->trace.keep_within_bounds)(E,in_coord_cc,in_coord_sph);
+ keep_within_bounds(E,in_coord_cc,in_coord_sph);
/* check whether tracer is within processor domain */
@@ -1097,10 +1255,10 @@
if (E->parallel.nprocz>1) icheck=icheck_processor_shell(E,j,in_coord_sph._rad);
if (icheck!=1) continue;
- if (E->parallel.nprocxy==1)
+ if (E->trace.full_tracers)
+ icheck=full_icheck_cap(E,0,in_coord_cc,in_coord_sph._rad);
+ else
icheck=regional_icheck_cap(E,0,in_coord_sph,in_coord_sph._rad);
- else
- icheck=full_icheck_cap(E,0,in_coord_cc,in_coord_sph._rad);
if (icheck==0) continue;
@@ -1237,7 +1395,7 @@
/* it is possible that if on phi=0 boundary, significant digits can push phi over 2pi */
- (E->trace.keep_within_bounds)(E,cc,sc);
+ keep_within_bounds(E,cc,sc);
Tracer new_tracer;
@@ -1505,7 +1663,7 @@
/********** ICHECK PROCESSOR SHELL *************/
-/* returns -99 if rad is below current shell */
+/* returns UNDEFINED_ELEMENT if rad is below current shell */
/* returns 0 if rad is above current shell */
/* returns 1 if rad is within current shell */
/* */
@@ -1529,7 +1687,7 @@
bottom_r = E->sx[j][3][1];
// First check bottom
- if (rad<bottom_r) return -99;
+ if (rad<bottom_r) return UNDEFINED_ELEMENT;
// Check top
if (rad<top_r) return 1;
@@ -1565,7 +1723,7 @@
/* nprocessor is right on bottom of me */
if (nprocessor == me-1) {
- if (icheck_processor_shell(E, j, rad) == -99) return 1;
+ if (icheck_processor_shell(E, j, rad) == UNDEFINED_ELEMENT) return 1;
else return 0;
}
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Viscosity_structures.c 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Viscosity_structures.c 2011-08-02 16:52:33 UTC (rev 18812)
@@ -1143,7 +1143,6 @@
void strain_rate_2_inv(struct All_variables *E, int m, float *EEDOT, int SQRT)
{
void get_rtf_at_ppts();
- void velo_from_element();
void construct_c3x3matrix_el();
void get_ba_p();
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/global_defs.h 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/global_defs.h 2011-08-02 16:52:33 UTC (rev 18812)
@@ -682,13 +682,6 @@
};
-struct CITCOM_GNOMONIC {
- /* gnomonic projected coordinate */
- double u;
- double v;
-};
-
-
#ifdef USE_HDF5
#include "hdf5_related.h"
#endif
@@ -726,9 +719,6 @@
/* for chemical convection & composition rheology */
struct COMPOSITION composition;
- struct CITCOM_GNOMONIC *gnomonic;
- double gnomonic_reference_phi;
-
struct COORD *eco[NCS];
struct IEN *ien[NCS]; /* global */
struct SIEN *sien[NCS];
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/prototypes.h
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/prototypes.h 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/prototypes.h 2011-08-02 16:52:33 UTC (rev 18812)
@@ -152,11 +152,10 @@
void make_regular_grid(struct All_variables *E);
void full_tracer_input(struct All_variables *);
void full_lost_souls(struct All_variables *);
-void full_get_shape_functions(struct All_variables *, double [9], int, SphericalCoord);
+void full_get_shape_functions(struct All_variables *, int &, double [6], int, SphericalCoord);
CartesianCoord full_get_velocity(struct All_variables *, int, int, SphericalCoord);
int full_icheck_cap(struct All_variables *, int, CartesianCoord, double);
int full_iget_element(struct All_variables *, int, int, CartesianCoord, SphericalCoord);
-void full_keep_within_bounds(struct All_variables *, CartesianCoord &, SphericalCoord &);
void analytical_test(struct All_variables *);
void analytical_runge_kutte(struct All_variables *, int, double, SphericalCoord &, CartesianCoord &, SphericalCoord &, CartesianCoord &, double *);
void analytical_test_function(struct All_variables *, SphericalCoord, SphericalCoord &, CartesianCoord &);
@@ -267,7 +266,7 @@
void assign_v_to_vector(struct All_variables *);
void v_from_vector_pseudo_surf(struct All_variables *);
void velo_from_element(struct All_variables *, float [4][9], int, int, int);
-void velo_from_element_d(struct All_variables *, double [4][9], int, int, int);
+void velo_from_element_d(struct All_variables *, CartesianCoord[9], int, int, int);
void p_to_nodes(struct All_variables *, double **, float **, int);
void visc_from_gint_to_nodes(struct All_variables *, float **, float **, int);
void visc_from_nodes_to_gint(struct All_variables *, float **, float **, int);
@@ -430,7 +429,6 @@
void regional_get_shape_functions(struct All_variables *, double [8], int, SphericalCoord);
double regional_interpolate_data(struct All_variables *, double [9], double [9]);
CartesianCoord regional_get_velocity(struct All_variables *, int, int, SphericalCoord);
-void regional_keep_within_bounds(struct All_variables *, CartesianCoord &, SphericalCoord &);
void regional_lost_souls(struct All_variables *);
/* Regional_version_dependent.c */
void regional_node_locations(struct All_variables *);
@@ -505,6 +503,7 @@
void count_tracers_of_flavors(struct All_variables *);
void initialize_tracers(struct All_variables *);
void get_neighboring_caps(struct All_variables *);
+void keep_within_bounds(struct All_variables *, CartesianCoord &, SphericalCoord &);
void allocate_tracer_arrays(struct All_variables *, int, int);
void expand_tracer_arrays(struct All_variables *, int);
void expand_later_array(struct All_variables *, int);
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h 2011-08-01 18:54:31 UTC (rev 18811)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h 2011-08-02 16:52:33 UTC (rev 18812)
@@ -75,10 +75,17 @@
return sqrt(xd*xd+yd*yd+zd*zd);
};
+ // Operations that return a new object
const CartesianCoord operator+(const CartesianCoord &other) const;
const CartesianCoord operator-(const CartesianCoord &other) const;
const CartesianCoord operator*(const double &val) const;
const CartesianCoord operator/(const double &val) const;
+
+ // Operations that modify this object
+ CartesianCoord & operator+=(const CartesianCoord &other);
+ CartesianCoord & operator-=(const CartesianCoord &other);
+ CartesianCoord & operator*=(const double &val);
+ CartesianCoord & operator/=(const double &val);
void operator=(const CartesianCoord &other) { _x = other._x; _y = other._y; _z = other._z; };
};
@@ -115,6 +122,14 @@
void setCartTrigBounds(unsigned int bnum, CartesianCoord cc, double cost, double sint, double cosf, double sinf);
};
+class TriElemLinearShapeFunc {
+private:
+ double sf[3][3];
+public:
+ TriElemLinearShapeFunc(CoordUV xy1, CoordUV xy2, CoordUV xy3);
+ double applyShapeFunc(const unsigned int i, const CoordUV uv) const;
+};
+
class Tracer {
private:
// Tracer position in spherical coordinates
@@ -227,6 +242,10 @@
/* Mesh information */
CapBoundary boundaries[13];
+ // Map of node numbers to gnomonic coordinates
+ std::map<int, CoordUV> *gnomonic;
+ double gnomonic_reference_phi;
+
/*********************/
/* for global model */
/*********************/
@@ -246,7 +265,7 @@
int *regtoel[13][5];
/* gnomonic shape functions */
- double *shape_coefs[13][3][10];
+ TriElemLinearShapeFunc **shape_coefs[13][2];
/**********************/
/* for regional model */
@@ -265,7 +284,4 @@
CartesianCoord (* get_velocity)(struct All_variables*, int, int,
SphericalCoord);
-
- void (* keep_within_bounds)(struct All_variables*,
- CartesianCoord &, SphericalCoord &);
};
More information about the CIG-COMMITS
mailing list