[cig-commits] r9271 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Feb 8 15:58:57 PST 2008
Author: tan2
Date: 2008-02-08 15:58:57 -0800 (Fri, 08 Feb 2008)
New Revision: 9271
Modified:
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
Log:
Refactored the functions.
* get shape functions on tracer location
* interpolate data using acquired shape functions
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2008-02-08 23:56:39 UTC (rev 9270)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2008-02-08 23:58:57 UTC (rev 9271)
@@ -843,10 +843,7 @@
return;
}
-
-/************************ GET VELOCITY ***************************************/
-/* */
-/* This function interpolates tracer velocity using gnominic interpolation. */
+/************************ GET SHAPE FUNCTION *********************************/
/* Real theta,phi,rad space is transformed into u,v space. This transformation */
/* maps great circles into straight lines. Here, elements boundaries are */
/* assumed to be great circle planes (not entirely true, it is actually only */
@@ -875,33 +872,28 @@
/* 5 6 5 7 */
/* 6 7 6 8 */
+void full_get_shape_functions(struct All_variables *E,
+ double shp[9], int nelem,
+ double theta, double phi, double rad)
+{
+ const int j = 1;
-void full_get_velocity(struct All_variables *E,
- int j, int nelem,
- double theta, double phi, double rad,
- double *velocity_vector)
-{
int iwedge,inum;
int i, kk;
int ival;
int itry;
- int sphere_key = 0;
double u,v;
double shape2d[4];
double shaperad[3];
double shape[7];
- double VV[4][9];
- double vx[7],vy[7],vz[7];
double x,y,z;
-
int maxlevel=E->mesh.levmax;
const double eps=-1e-4;
void sphere_to_cart();
- void velo_from_element_d();
/* find u and v using spherical coordinates */
@@ -982,18 +974,76 @@
/* Sum of shape functions is 1 */
- shape[1]=shaperad[1]*shape2d[1];
- shape[2]=shaperad[1]*shape2d[2];
- shape[3]=shaperad[1]*shape2d[3];
- shape[4]=shaperad[2]*shape2d[1];
- shape[5]=shaperad[2]*shape2d[2];
- shape[6]=shaperad[2]*shape2d[3];
+ 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];
/** debug **
fprintf(E->trace.fpt, "shp: %e %e %e %e %e %e\n",
- shape[1], shape[2], shape[3], shape[4], shape[5], shape[6]);
+ shp[1], shp[2], shp[3], shp[4], shp[5], shp[6]);
/**/
+ return;
+}
+
+
+double full_interpolate_data(struct All_variables *E,
+ double shp[9], double data[9])
+{
+ int iwedge = shp[0];
+
+ if (iwedge==1)
+ return data[1]*shp[1] + data[2]*shp[2] + shp[3]*data[3]
+ + data[6]*shp[4] + data[6]*shp[5] + shp[7]*data[6];
+
+ if (iwedge==2)
+ return data[1]*shp[1] + data[3]*shp[2] + shp[4]*data[3]
+ + data[5]*shp[4] + data[7]*shp[5] + shp[8]*data[6];
+}
+
+
+/************************ GET VELOCITY ***************************************/
+/* */
+/* This function interpolates tracer velocity using gnominic interpolation. */
+/* The element is divided into 2 wedges in which standard shape functions */
+/* are used to interpolate velocity. */
+/* */
+/* Wedge information: */
+/* */
+/* Wedge 1 Wedge 2 */
+/* _______ _______ */
+/* */
+/* wedge_node real_node wedge_node real_node */
+/* ---------- --------- ---------- --------- */
+/* */
+/* 1 1 1 1 */
+/* 2 2 2 3 */
+/* 3 3 3 4 */
+/* 4 5 4 5 */
+/* 5 6 5 7 */
+/* 6 7 6 8 */
+
+void full_get_velocity(struct All_variables *E,
+ int j, int nelem,
+ double theta, double phi, double rad,
+ double *velocity_vector)
+{
+ int iwedge;
+ const int sphere_key = 0;
+
+ double shape[9];
+ double VV[4][9];
+ double vx[7],vy[7],vz[7];
+
+ void velo_from_element_d();
+
+ full_get_shape_functions(E, shape, nelem, theta, phi, rad);
+ iwedge=shape[0];
+
/* get cartesian velocity */
velo_from_element_d(E, VV, j, nelem, sphere_key);
Modified: mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c 2008-02-08 23:56:39 UTC (rev 9270)
+++ mc/3D/CitcomS/trunk/lib/Regional_tracer_advection.c 2008-02-08 23:58:57 UTC (rev 9271)
@@ -421,9 +421,9 @@
}
-static void get_shape_functions(struct All_variables *E,
- double w[9], int nelem,
- double theta, double phi, double rad)
+void regional_get_shape_functions(struct All_variables *E,
+ double shp[9], int nelem,
+ double theta, double phi, double rad)
{
int e, i, j, k;
int elx, ely, elz;
@@ -473,29 +473,41 @@
/* compute volumetic weighting functions */
volume = dx*dz*dy;
- w[1] = (dx-tr_dx) * (dy-tr_dy) * (dz-tr_dz) / volume;
- w[2] = tr_dx * (dy-tr_dy) * (dz-tr_dz) / volume;
- w[3] = tr_dx * tr_dy * (dz-tr_dz) / volume;
- w[4] = (dx-tr_dx) * tr_dy * (dz-tr_dz) / volume;
- w[5] = (dx-tr_dx) * (dy-tr_dy) * tr_dz / volume;
- w[6] = tr_dx * (dy-tr_dy) * tr_dz / volume;
- w[7] = tr_dx * tr_dy * tr_dz / volume;
- w[8] = (dx-tr_dx) * tr_dy * tr_dz / volume;
+ shp[1] = (dx-tr_dx) * (dy-tr_dy) * (dz-tr_dz) / volume;
+ shp[2] = tr_dx * (dy-tr_dy) * (dz-tr_dz) / volume;
+ shp[3] = tr_dx * tr_dy * (dz-tr_dz) / volume;
+ shp[4] = (dx-tr_dx) * tr_dy * (dz-tr_dz) / volume;
+ shp[5] = (dx-tr_dx) * (dy-tr_dy) * tr_dz / volume;
+ shp[6] = tr_dx * (dy-tr_dy) * tr_dz / volume;
+ shp[7] = tr_dx * tr_dy * tr_dz / volume;
+ shp[8] = (dx-tr_dx) * tr_dy * tr_dz / volume;
/** debug **
fprintf(E->trace.fpt, "dr=(%e,%e,%e) tr_dr=(%e,%e,%e)\n",
dx, dy, dz, tr_dx, tr_dy, tr_dz);
fprintf(E->trace.fpt, "shp: %e %e %e %e %e %e %e %e\n",
- w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8]);
+ shp[1], shp[2], shp[3], shp[4], shp[5], shp[6], shp[7], shp[8]);
fprintf(E->trace.fpt, "sum(shp): %e\n",
- w[1]+ w[2]+ w[3]+ w[4]+ w[5]+ w[6]+ w[7]+ w[8]);
+ shp[1]+ shp[2]+ shp[3]+ shp[4]+ shp[5]+ shp[6]+ shp[7]+ shp[8]);
fflush(E->trace.fpt);
/**/
return;
}
+double regional_interpolate_data(struct All_variables *E,
+ double shp[9], double data[9])
+{
+ int n;
+ double result = 0;
+ for(n=1; n<=8; n++)
+ result += data[n] * shp[n];
+
+ return result;
+}
+
+
/******** GET VELOCITY ***************************************/
void regional_get_velocity(struct All_variables *E,
@@ -505,12 +517,12 @@
{
void velo_from_element_d();
- double weight[9], VV[4][9], tmp;
+ double shp[9], VV[4][9], tmp;
int n, d, node;
const int sphere_key = 0;
/* get shape functions at (theta, phi, rad) */
- get_shape_functions(E, weight, nelem, theta, phi, rad);
+ regional_get_shape_functions(E, shp, nelem, theta, phi, rad);
/* get cartesian velocity */
@@ -527,7 +539,7 @@
for(d=1; d<=3; d++) {
for(n=1; n<=8; n++)
- velocity_vector[d] += VV[d][n] * weight[n];
+ velocity_vector[d] += VV[d][n] * shp[n];
}
@@ -541,7 +553,7 @@
tmp = 0;
for(n=1; n<=8; n++)
- tmp += E->sx[m][1][E->ien[m][nelem].node[n]] * weight[n];
+ tmp += E->sx[m][1][E->ien[m][nelem].node[n]] * shp[n];
fprintf(E->trace.fpt, "THETA: %e -> %e\n", theta, tmp);
More information about the cig-commits
mailing list