[cig-commits] r9110 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Tue Jan 22 11:38:43 PST 2008
Author: tan2
Date: 2008-01-22 11:38:43 -0800 (Tue, 22 Jan 2008)
New Revision: 9110
Modified:
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Obsolete.c
mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
mc/3D/CitcomS/trunk/lib/Size_does_matter.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Computes the derivatives of shape functions and jacobians only once and stores them for later use. This reduces the total cpu time of cookbook8 by 13%.
- get_global_shape_fn_sph() computes and stores the derivatives of shape functions and jacobians.
- Moved get_global_shape_fn() to Obsolete.c
- Calling construct_shape_function_derivatives() in initial_mesh_solver_setup(). This function call get_global_shape_fn_sph() for each element.
- get_rtf_at_vpts() and get_rtf_at_ppts() are for coord. transformation matrix at vpts and ppts respectively
- Passed arguments by references in PG solvers instead of by values to avoid copying.
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -44,7 +44,7 @@
double **fielddot, double **Dfielddot);
static void pg_solver(struct All_variables *E,
double **T, double **Tdot, double **DTdot,
- struct SOURCES Q0,
+ struct SOURCES *Q0,
double diff, int bc, unsigned int **FLAGS);
static void pg_shape_fn(struct All_variables *E, int el,
struct Shape_function *PG,
@@ -52,12 +52,12 @@
float VV[4][9], double rtf[4][9],
double diffusion, int m);
static void element_residual(struct All_variables *E, int el,
- struct Shape_function PG,
- struct Shape_function_dx GNx,
- struct Shape_function_dA dOmega,
+ struct Shape_function *PG,
+ struct Shape_function_dx *GNx,
+ struct Shape_function_dA *dOmega,
float VV[4][9],
double **field, double **fielddot,
- struct SOURCES Q0,
+ struct SOURCES *Q0,
double Eres[9], double rtf[4][9],
double diff, float **BC,
unsigned int **FLAGS, int m);
@@ -246,7 +246,7 @@
process_heating(E, psc_pass);
/* XXX: replace inputdiff with refstate.thermal_conductivity */
- pg_solver(E,E->T,E->Tdot,DTdot,E->convection.heat_sources,E->control.inputdiff,1,E->node);
+ pg_solver(E,E->T,E->Tdot,DTdot,&(E->convection.heat_sources),E->control.inputdiff,1,E->node);
corrector(E,E->T,E->Tdot,DTdot);
temperatures_conform_bcs(E);
}
@@ -389,10 +389,10 @@
static void pg_solver(struct All_variables *E,
double **T, double **Tdot, double **DTdot,
- struct SOURCES Q0,
+ struct SOURCES *Q0,
double diff, int bc, unsigned int **FLAGS)
{
- void get_global_shape_fn();
+ void get_rtf_vpts();
void velo_from_element();
int el,e,a,i,a1,m;
@@ -400,14 +400,12 @@
float VV[4][9];
struct Shape_function PG;
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
const int dims=E->mesh.nsd;
const int dofs=E->mesh.dof;
const int ends=enodes[dims];
const int sphere_key = 1;
+ const int lev=E->mesh.levmax;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++)
@@ -418,13 +416,13 @@
velo_from_element(E,VV,m,el,sphere_key);
- get_global_shape_fn(E, el, &GN, &GNx, &dOmega, 0,
- sphere_key, rtf, E->mesh.levmax, m);
+ get_rtf_vpts(E, m, lev, el, rtf);
/* XXX: replace diff with refstate.thermal_conductivity */
- pg_shape_fn(E, el, &PG, &GNx, VV,
+ pg_shape_fn(E, el, &PG, &(E->gNX[m][el]), VV,
rtf, diff, m);
- element_residual(E, el, PG, GNx, dOmega, VV, T, Tdot,
+ element_residual(E, el, &PG, &(E->gNX[m][el]), &(E->gDA[m][el]),
+ VV, T, Tdot,
Q0, Eres, rtf, diff, E->sphere.cap[m].TB,
FLAGS, m);
@@ -435,7 +433,7 @@
} /* next element */
- (E->exchange_node_d)(E,DTdot,E->mesh.levmax);
+ (E->exchange_node_d)(E,DTdot,lev);
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++) {
@@ -525,12 +523,12 @@
========================================= */
static void element_residual(struct All_variables *E, int el,
- struct Shape_function PG,
- struct Shape_function_dx GNx,
- struct Shape_function_dA dOmega,
+ struct Shape_function *PG,
+ struct Shape_function_dx *GNx,
+ struct Shape_function_dA *dOmega,
float VV[4][9],
double **field, double **fielddot,
- struct SOURCES Q0,
+ struct SOURCES *Q0,
double Eres[9], double rtf[4][9],
double diff, float **BC,
unsigned int **FLAGS, int m)
@@ -579,9 +577,9 @@
for(i=1;i<=vpts;i++) {
dT[i] += DT * E->N.vpt[GNVINDEX(j,i)];
- tx1[i] += GNx.vpt[GNVXINDEX(0,j,i)] * T * rtf[3][i];
- tx2[i] += GNx.vpt[GNVXINDEX(1,j,i)] * T * sint[i];
- tx3[i] += GNx.vpt[GNVXINDEX(2,j,i)] * T;
+ tx1[i] += GNx->vpt[GNVXINDEX(0,j,i)] * T * rtf[3][i];
+ tx2[i] += GNx->vpt[GNVXINDEX(1,j,i)] * T * sint[i];
+ tx3[i] += GNx->vpt[GNVXINDEX(2,j,i)] * T;
sfn = E->N.vpt[GNVINDEX(j,i)];
v1[i] += VV[1][j] * sfn;
v2[i] += VV[2][j] * sfn;
@@ -591,7 +589,7 @@
/* Q=0.0;
for(i=0;i<Q0.number;i++)
- Q += Q0.Q[i] * exp(-Q0.lambda[i] * (E->monitor.elapsed_time+Q0.t_offset));
+ Q += Q0->Q[i] * exp(-Q0->lambda[i] * (E->monitor.elapsed_time+Q0->t_offset));
*/
/* heat production */
@@ -626,13 +624,13 @@
Eres[j]=0.0;
for(i=1;i<=vpts;i++)
Eres[j] -=
- PG.vpt[GNVINDEX(j,i)] * dOmega.vpt[i]
+ PG->vpt[GNVINDEX(j,i)] * dOmega->vpt[i]
* ((dT[i] + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i])*rho*cp
- heating )
- + diff * dOmega.vpt[i] * E->heating_latent[m][el]
- * (GNx.vpt[GNVXINDEX(0,j,i)]*tx1[i]*rtf[3][i] +
- GNx.vpt[GNVXINDEX(1,j,i)]*tx2[i]*sint[i] +
- GNx.vpt[GNVXINDEX(2,j,i)]*tx3[i] );
+ + diff * dOmega->vpt[i] * E->heating_latent[m][el]
+ * (GNx->vpt[GNVXINDEX(0,j,i)]*tx1[i]*rtf[3][i] +
+ GNx->vpt[GNVXINDEX(1,j,i)]*tx2[i]*sint[i] +
+ GNx->vpt[GNVXINDEX(2,j,i)]*tx3[i] );
}
}
@@ -640,7 +638,7 @@
for(j=1;j<=ends;j++) {
Eres[j]=0.0;
for(i=1;i<=vpts;i++)
- Eres[j] -= PG.vpt[GNVINDEX(j,i)] * dOmega.vpt[i]
+ Eres[j] -= PG->vpt[GNVINDEX(j,i)] * dOmega->vpt[i]
* (dT[i] - heating + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i]);
}
}
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -311,11 +311,8 @@
const double two = 2.0;
const double two_thirds = 2.0/3.0;
- void get_global_shape_fn();
+ void get_rtf_vpts();
void construct_c3x3matrix_el();
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
double ba[9][9][4][7]; /* integration points,node,3x6 matrix */
@@ -323,9 +320,8 @@
const int vpts = VPOINTS3D;
const int ends = ENODES3D;
const int dims=E->mesh.nsd;
- const int sphere_key = 1;
- get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,lev,m);
+ get_rtf_vpts(E, m, lev, el, rtf);
if (iconv || (el-1)%E->lmesh.ELZ[lev]==0)
construct_c3x3matrix_el(E,el,&E->element_Cc,&E->element_Ccx,lev,m,0);
@@ -334,10 +330,10 @@
quadrature point. Nx[d] is the derivative wrt x[d]. */
for(k=1;k<=vpts;k++) {
- W[k]=g_point[k].weight[dims-1]*dOmega.vpt[k]*E->EVI[lev][m][(el-1)*vpts+k];
+ W[k]=g_point[k].weight[dims-1]*E->GDA[lev][m][el].vpt[k]*E->EVI[lev][m][(el-1)*vpts+k];
}
- get_ba(&(E->N), &GNx, &E->element_Cc, &E->element_Ccx,
+ get_ba(&(E->N), &(E->GNX[lev][m][el]), &E->element_Cc, &E->element_Ccx,
rtf, E->mesh.nsd, ba);
for(a=1;a<=ends;a++)
@@ -826,26 +822,19 @@
void get_elt_c(struct All_variables *E, int el,
higher_precision elt_c[24][1], int lev, int m)
{
- void get_global_shape_fn();
void construct_c3x3matrix_el();
int p, a, i, j, nz;
double temp, beta, rho_avg, x[4];
- struct Shape_function GN;
- struct Shape_function_dx GNx;
- struct Shape_function_dA dOmega;
- double rtf[4][9], rho[9];
+ double rho[9];
const int dims = E->mesh.nsd;
const int ends = enodes[dims];
- const int sphere_key = 1;
- get_global_shape_fn(E,el,&GN,&GNx,&dOmega,2,sphere_key,rtf,lev,m);
-
if ((el-1)%E->lmesh.ELZ[lev]==0)
construct_c3x3matrix_el(E,el,&E->element_Cc,&E->element_Ccx,lev,m,1);
- temp = p_point[1].weight[dims-1] * dOmega.ppt[1];
+ temp = p_point[1].weight[dims-1] * E->GDA[lev][m][el].ppt[1];
switch (E->refstate.choice) {
case 1:
@@ -882,7 +871,7 @@
for(a=1;a<=ends;a++) {
for (i=1;i<=dims;i++) {
- x[i] = rho[a] * GNx.ppt[GNPXINDEX(2,a,1)]
+ x[i] = rho[a] * E->GNX[lev][m][el].ppt[GNPXINDEX(2,a,1)]
* E->N.ppt[GNPINDEX(a,1)]
* E->element_Cc.ppt[BPINDEX(3,i,a,1)];
}
@@ -908,44 +897,37 @@
higher_precision elt_del[24][1];
int lev;
-{ void get_global_shape_fn();
+{
+ void get_rtf_ppts();
void construct_c3x3matrix_el();
int p,a,i;
double ra,ct,si,x[4],rtf[4][9];
double temp;
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
-
const int dims=E->mesh.nsd;
const int ends=enodes[dims];
- const int sphere_key = 1;
/* Special case, 4/8 node bilinear cartesian square/cube element -> 1 pressure point */
-/* es = (el-1)/E->lmesh.ELZ[lev]+1; */
-
if ((el-1)%E->lmesh.ELZ[lev]==0)
construct_c3x3matrix_el(E,el,&E->element_Cc,&E->element_Ccx,lev,m,1);
- get_global_shape_fn(E,el,&GN,&GNx,&dOmega,2,sphere_key,rtf,lev,m);
+ get_rtf_ppts(E, m, lev, el, rtf);
+ temp=p_point[1].weight[dims-1] * E->GDA[lev][m][el].ppt[1];
- temp=p_point[1].weight[dims-1] * dOmega.ppt[1];
-
ra = rtf[3][1];
si = 1.0/sin(rtf[1][1]);
ct = cos(rtf[1][1])*si;
for(a=1;a<=ends;a++) {
for (i=1;i<=dims;i++)
- x[i]=GNx.ppt[GNPXINDEX(2,a,1)]*E->element_Cc.ppt[BPINDEX(3,i,a,1)]
+ x[i]=E->GNX[lev][m][el].ppt[GNPXINDEX(2,a,1)]*E->element_Cc.ppt[BPINDEX(3,i,a,1)]
+ 2.0*ra*E->N.ppt[GNPINDEX(a,1)]*E->element_Cc.ppt[BPINDEX(3,i,a,1)]
- + ra*(GNx.ppt[GNPXINDEX(0,a,1)]*E->element_Cc.ppt[BPINDEX(1,i,a,1)]
+ + ra*(E->GNX[lev][m][el].ppt[GNPXINDEX(0,a,1)]*E->element_Cc.ppt[BPINDEX(1,i,a,1)]
+E->N.ppt[GNPINDEX(a,1)]*E->element_Ccx.ppt[BPXINDEX(1,i,1,a,1)]
+ct*E->N.ppt[GNPINDEX(a,1)]*E->element_Cc.ppt[BPINDEX(1,i,a,1)]
- +si*(GNx.ppt[GNPXINDEX(1,a,1)]*E->element_Cc.ppt[BPINDEX(2,i,a,1)]
+ +si*(E->GNX[lev][m][el].ppt[GNPXINDEX(1,a,1)]*E->element_Cc.ppt[BPINDEX(2,i,a,1)]
+E->N.ppt[GNPINDEX(a,1)]*E->element_Ccx.ppt[BPXINDEX(2,i,2,a,1)]));
p=dims*(a-1);
@@ -953,7 +935,7 @@
elt_del[p+1][0] = -x[2] * temp;
elt_del[p+2][0] = -x[3] * temp;
- /* fprintf (E->fp,"B= %d %d %g %g %g %g %g\n",el,a,dOmega.ppt[1],GNx.ppt[GNPXINDEX(0,a,1)],GNx.ppt[GNPXINDEX(1,a,1)],elt_del[p][0],elt_del[p+1][0]);
+ /* fprintf (E->fp,"B= %d %d %g %g %g %g %g\n",el,a,E->GDA[lev][m][el].ppt[1],E->GNX[lev][m][el].ppt[GNPXINDEX(0,a,1)],E->GNX[lev][m][el].ppt[GNPXINDEX(1,a,1)],elt_del[p][0],elt_del[p+1][0]);
*/
}
@@ -978,25 +960,16 @@
const unsigned int vbc_flag[] = {0, VBX, VBY, VBZ};
double force[9],force_at_gs[9],elt_k[24*24];
- double rtf[4][9];
- void get_global_shape_fn();
void construct_c3x3matrix_el();
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
-
const int dims=E->mesh.nsd;
const int n=loc_mat_size[dims];
const int ends=enodes[dims];
const int vpts=vpoints[dims];
- const int sphere_key=1;
es = (el-1)/E->lmesh.elz + 1;
- get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,E->mesh.levmax,m);
-
if ((el-1)%E->lmesh.elz==0)
construct_c3x3matrix_el(E,el,&E->element_Cc,&E->element_Ccx,E->mesh.levmax,m,0);
@@ -1018,7 +991,7 @@
for(j=1;j<=vpts;j++) /*compute sum(Na(j)*F(j)*det(j)) */
elt_f[p] += force_at_gs[j] * E->N.vpt[GNVINDEX(a,j)]
- *dOmega.vpt[j]*g_point[j].weight[dims-1]
+ *E->gDA[m][el].vpt[j]*g_point[j].weight[dims-1]
*E->element_Cc.vpt[BVINDEX(3,i,a,j)];
/* imposed velocity terms */
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -312,21 +312,11 @@
int average;
{
- void get_global_shape_fn();
- void float_global_operation();
-
- double rtf[4][9];
-
int n,i,j,k,el,m;
float volume,integral,volume1,integral1;
- struct Shape_function GN;
- struct Shape_function_dx GNx;
- struct Shape_function_dA dOmega;
-
const int vpts = vpoints[E->mesh.nsd];
const int ends = enodes[E->mesh.nsd];
- const int sphere_key=1;
volume1=0.0;
integral1=0.0;
@@ -334,13 +324,11 @@
for (m=1;m<=E->sphere.caps_per_proc;m++)
for (el=1;el<=E->lmesh.nel;el++) {
- get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,E->mesh.levmax,m);
-
for(j=1;j<=vpts;j++)
for(i=1;i<=ends;i++) {
n = E->ien[m][el].node[i];
- volume1 += E->N.vpt[GNVINDEX(i,j)] * dOmega.vpt[j];
- integral1 += Z[m][n] * E->N.vpt[GNVINDEX(i,j)] * dOmega.vpt[j];
+ volume1 += E->N.vpt[GNVINDEX(i,j)] * E->gDA[m][el].vpt[j];
+ integral1 += Z[m][n] * E->N.vpt[GNVINDEX(i,j)] * E->gDA[m][el].vpt[j];
}
}
@@ -368,19 +356,11 @@
int average;
{
- void get_global_shape_fn();
-
- double rtf[4][9];
int n,i,j,el,m;
double volume,integral,volume1,integral1;
- struct Shape_function GN;
- struct Shape_function_dx GNx;
- struct Shape_function_dA dOmega;
-
const int vpts = vpoints[E->mesh.nsd];
const int ends = enodes[E->mesh.nsd];
- const int sphere_key=1;
volume1=0.0;
integral1=0.0;
@@ -388,12 +368,11 @@
for (m=1;m<=E->sphere.caps_per_proc;m++)
for (el=1;el<=E->lmesh.nel;el++) {
- get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,E->mesh.levmax,m);
for(j=1;j<=vpts;j++)
for(i=1;i<=ends;i++) {
n = E->ien[m][el].node[i];
- volume1 += E->N.vpt[GNVINDEX(i,j)] * dOmega.vpt[j];
- integral1 += Z[m][n] * E->N.vpt[GNVINDEX(i,j)] * dOmega.vpt[j];
+ volume1 += E->N.vpt[GNVINDEX(i,j)] * E->gDA[m][el].vpt[j];
+ integral1 += Z[m][n] * E->N.vpt[GNVINDEX(i,j)] * E->gDA[m][el].vpt[j];
}
}
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -57,6 +57,7 @@
void check_bc_consistency(struct All_variables*);
void construct_elt_gs(struct All_variables*);
void construct_elt_cs(struct All_variables*);
+void construct_shape_function_derivatives(struct All_variables *E);
void construct_id(struct All_variables*);
void construct_ien(struct All_variables*);
void construct_lm(struct All_variables*);
@@ -141,6 +142,7 @@
construct_sub_element(E);
construct_shape_functions(E);
+ construct_shape_function_derivatives(E);
construct_elt_gs(E);
if(E->control.inv_gruneisen != 0)
construct_elt_cs(E);
@@ -644,6 +646,9 @@
for (k=1;k<=4;k++)
E->sphere.angle1[i][j][k] = (double *) malloc((snel+1)*sizeof(double));
+ E->GNX[i][j] = (struct Shape_function_dx *)malloc((nel+1)*sizeof(struct Shape_function_dx));
+ E->GDA[i][j] = (struct Shape_function_dA *)malloc((nel+1)*sizeof(struct Shape_function_dA));
+
E->MASS[i][j] = (float *) malloc((nno+1)*sizeof(float));
E->ECO[i][j] = (struct COORD *) malloc((nno+2)*sizeof(struct COORD));
@@ -979,6 +984,8 @@
E->ccx[j] = E->CCX[E->mesh.levmax][j];
E->Mass[j] = E->MASS[E->mesh.levmax][j];
E->element[j] = E->ELEMENT[E->mesh.levmax][j];
+ E->gDA[j] = E->GDA[E->mesh.levmax][j];
+ E->gNX[j] = E->GNX[E->mesh.levmax][j];
for (i=1;i<=E->mesh.nsd;i++) {
E->x[j][i] = E->X[E->mesh.levmax][j][i];
Modified: mc/3D/CitcomS/trunk/lib/Obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Obsolete.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Obsolete.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -38,6 +38,162 @@
/* =========================================================== */
+
+/* ==================================================================================
+ Function to give the global shape function from the local: Assumes ORTHOGONAL MESH
+ ================================================================================== */
+
+void get_global_shape_fn(E,el,GN,GNx,dOmega,pressure,sphere,rtf,lev,m)
+ struct All_variables *E;
+ int el,m;
+ struct Shape_function *GN;
+ struct Shape_function_dx *GNx;
+ struct Shape_function_dA *dOmega;
+ int pressure,lev,sphere;
+ double rtf[4][9];
+{
+ int i,j,k,d,e;
+ double jacobian;
+ double determinant();
+ double cofactor(),myatan();
+ void form_rtf_bc();
+
+ struct Shape_function_dx LGNx;
+
+ double dxda[4][4],cof[4][4],x[4],bc[4][4];
+
+
+ const int dims=E->mesh.nsd;
+ const int ends=enodes[dims];
+ const int vpts=vpoints[dims];
+ const int ppts=ppoints[dims];
+
+
+ if(pressure < 2) {
+ 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);
+ dOmega->vpt[k] = jacobian;
+
+ for(d=1;d<=dims;d++)
+ for(e=1;e<=dims;e++)
+ cof[d][e]=cofactor(dxda,d,e,dims);
+
+ if (sphere) {
+
+ 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++) {
+ GNx->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)];
+ }
+ }
+ else {
+ for(j=1;j<=ends;j++)
+ for(d=1;d<=dims;d++) {
+ GNx->vpt[GNVXINDEX(d-1,j,k)] = 0.0;
+ for(e=1;e<=dims;e++)
+ GNx->vpt[GNVXINDEX(d-1,j,k)] +=
+ E->Nx.vpt[GNVXINDEX(e-1,j,k)] *cof[e][d];
+
+ GNx->vpt[GNVXINDEX(d-1,j,k)] /= jacobian;
+ }
+ }
+ } /* end for k */
+ } /* end for pressure */
+
+ if(pressure > 0 && pressure < 3) {
+ 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);
+ dOmega->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);
+
+ if (sphere) {
+ 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++) {
+ GNx->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)];
+ }
+ }
+
+ else {
+ for(j=1;j<=ends;j++)
+ for(d=1;d<=dims;d++) {
+ GNx->ppt[GNPXINDEX(d-1,j,k)]=0.0;
+ for(e=1;e<=dims;e++)
+ GNx->ppt[GNPXINDEX(d-1,j,k)] +=
+ E->Nx.ppt[GNPXINDEX(e-1,j,k)]*cof[e][d];
+ GNx->ppt[GNPXINDEX(d-1,j,k)] /= jacobian;
+ }
+ }
+
+ } /* end for k int */
+ } /* end for pressure */
+
+
+ return;
+}
+
+
void get_global_1d_shape_fn_1(E,el,GM,dGammax,nodal,m)
struct All_variables *E;
int el,nodal,m;
Modified: mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -56,12 +56,7 @@
float *flux[NCS],*SU[NCS],*RU[NCS];
float VV[4][9],u[9],T[9],dTdz[9],area,uT;
float *sum_h;
- double rtf[4][9];
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
- void get_global_shape_fn();
void velo_from_element();
void sum_across_surface();
void return_horiz_ave();
@@ -75,7 +70,6 @@
const int lev = E->mesh.levmax;
const int sphere_key=1;
-
sum_h = (float *) malloc((5)*sizeof(float));
for(i=0;i<=4;i++)
sum_h[i] = 0.0;
@@ -89,7 +83,6 @@
}
for(e=1;e<=E->lmesh.nel;e++) {
- get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,sphere_key,rtf,lev,m);
velo_from_element(E,VV,m,e,sphere_key);
@@ -100,7 +93,7 @@
for(j=1;j<=ends;j++) {
u[i] += VV[3][j]*E->N.vpt[GNVINDEX(j,i)];
T[i] += E->T[m][E->ien[m][e].node[j]]*E->N.vpt[GNVINDEX(j,i)];
- dTdz[i] += -E->T[m][E->ien[m][e].node[j]]*GNx.vpt[GNVXINDEX(2,j,i)];
+ dTdz[i] += -E->T[m][E->ien[m][e].node[j]]*E->gNX[m][e].vpt[GNVXINDEX(2,j,i)];
}
}
@@ -108,7 +101,7 @@
area = 0.0;
for(i=1;i<=vpts;i++) {
/* XXX: missing unit conversion, heat capacity and thermal conductivity */
- uT += u[i]*T[i]*dOmega.vpt[i] + dTdz[i]*dOmega.vpt[i];
+ uT += u[i]*T[i]*E->gDA[m][e].vpt[i] + dTdz[i]*E->gDA[m][e].vpt[i];
}
uT /= E->eco[m][e].area;
Modified: mc/3D/CitcomS/trunk/lib/Size_does_matter.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Size_does_matter.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Size_does_matter.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -42,189 +42,219 @@
return; }
+/* ======================================================================
+ ====================================================================== */
-/* ==================================================================================
- Function to give the global shape function from the local: Assumes ORTHOGONAL MESH
- ================================================================================== */
+static void form_rtf_bc(int k, double x[4],
+ double rtf[4][9], double bc[4][4])
+{
+ double myatan();
-void get_global_shape_fn(E,el,GN,GNx,dOmega,pressure,sphere,rtf,lev,m)
- struct All_variables *E;
- int el,m;
- struct Shape_function *GN;
- struct Shape_function_dx *GNx;
- struct Shape_function_dA *dOmega;
- int pressure,lev,sphere;
- double rtf[4][9];
+ 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]);
+
+ bc[1][1] = x[3]*cos(rtf[2][k]);
+ bc[1][2] = x[3]*sin(rtf[2][k]);
+ bc[1][3] = -sin(rtf[1][k])/rtf[3][k];
+ bc[2][1] = -x[2];
+ 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;
- double determinant();
- double cofactor(),myatan();
- void form_rtf_bc();
+ int i,j,k,d,e;
+ double jacobian;
+ double determinant();
+ double cofactor(),myatan();
+ void form_rtf_bc();
- struct Shape_function_dx LGNx;
+ struct Shape_function_dx LGNx;
- double dxda[4][4],cof[4][4],x[4],bc[4][4];
+ 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;
- const int dims=E->mesh.nsd;
- const int ends=enodes[dims];
- const int vpts=vpoints[dims];
- const int ppts=ppoints[dims];
-
- if(pressure < 2) {
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++) {
+ 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(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)];
+ 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);
- dOmega->vpt[k] = jacobian;
+ 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);
+ for(d=1;d<=dims;d++)
+ for(e=1;e<=dims;e++)
+ cof[d][e]=cofactor(dxda,d,e,dims);
- if (sphere) {
-
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];
+ 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;
+ LGNx.vpt[GNVXINDEX(d-1,j,k)] /= jacobian;
}
for(j=1;j<=ends;j++)
- for(d=1;d<=dims;d++) {
- GNx->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)];
+ 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)];
}
- }
- else {
- for(j=1;j<=ends;j++)
- for(d=1;d<=dims;d++) {
- GNx->vpt[GNVXINDEX(d-1,j,k)] = 0.0;
+ } /* 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++)
- GNx->vpt[GNVXINDEX(d-1,j,k)] +=
- E->Nx.vpt[GNVXINDEX(e-1,j,k)] *cof[e][d];
-
- GNx->vpt[GNVXINDEX(d-1,j,k)] /= jacobian;
- }
+ dxda[d][e]=0.0;
}
- } /* end for k */
- } /* end for pressure */
- if(pressure > 0 && pressure < 3) {
- 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(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)];
- 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;
- jacobian = determinant(dxda,E->mesh.nsd);
- dOmega->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);
- for(d=1;d<=dims;d++)
- for(e=1;e<=dims;e++)
- cof[d][e]=cofactor(dxda,d,e,E->mesh.nsd);
-
- if (sphere) {
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(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++) {
- GNx->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)];
- }
- }
+ 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)];
+ }
- else {
- for(j=1;j<=ends;j++)
- for(d=1;d<=dims;d++) {
- GNx->ppt[GNPXINDEX(d-1,j,k)]=0.0;
- for(e=1;e<=dims;e++)
- GNx->ppt[GNPXINDEX(d-1,j,k)] +=
- E->Nx.ppt[GNPXINDEX(e-1,j,k)]*cof[e][d];
- GNx->ppt[GNPXINDEX(d-1,j,k)] /= jacobian;
+ } /* 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);
}
- }
- } /* end for k int */
- } /* end for pressure */
+ return;
+}
- return;
+void get_rtf_vpts(struct All_variables *E, int m, int lev, int el,
+ double rtf[4][9])
+{
+ int i, k, d;
+ double x[4];
+
+ double myatan();
+
+ 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]);
+ rtf[1][k] = acos(x[3]*rtf[3][k]);
+ rtf[2][k] = myatan(x[2],x[1]);
+ }
+
+ return;
}
-/* ======================================================================
- ====================================================================== */
-void form_rtf_bc(k,x,rtf,bc)
- int k;
- double x[4],rtf[4][9],bc[4][4];
- {
+void get_rtf_ppts(struct All_variables *E, int m, int lev, int el,
+ double rtf[4][9])
+{
+ int i, k, d;
+ double x[4];
- double myatan();
+ double myatan();
- 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]);
+ const int dims = E->mesh.nsd;
+ const int ends = ENODES3D;
+ const int ppts = PPOINTS3D;
- bc[1][1] = x[3]*cos(rtf[2][k]);
- bc[1][2] = x[3]*sin(rtf[2][k]);
- bc[1][3] = -sin(rtf[1][k])/rtf[3][k];
- bc[2][1] = -x[2];
- 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];
+ for(k=1;k<=ppts;k++) { /* all of the ppoints */
+ for(d=1;d<=dims;d++)
+ x[d]=0.0;
- return;
- }
+ 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)
{
@@ -871,14 +901,9 @@
{
int m,node,i,nint,e,lev;
int n[9], nz;
- void get_global_shape_fn();
- double myatan(),rtf[4][9],area,centre[4],temp[9],temp2[9],dx1,dx2,dx3;
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
+ double myatan(),area,centre[4],temp[9],temp2[9],dx1,dx2,dx3;
const int vpts=vpoints[E->mesh.nsd];
- const int sphere_key=1;
/* ECO .size can also be defined here */
@@ -890,8 +915,6 @@
for(e=1;e<=E->lmesh.NEL[lev];e++) {
- get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,sphere_key,rtf,lev,m);
-
area = centre[1] = centre[2] = centre[3] = 0.0;
for(node=1;node<=enodes[E->mesh.nsd];node++)
@@ -949,13 +972,13 @@
/* volume (area in 2D) of this element */
for(nint=1;nint<=vpts;nint++)
- area += g_point[nint].weight[E->mesh.nsd-1] * dOmega.vpt[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] += dOmega.vpt[nint]*g_point[nint].weight[E->mesh.nsd-1]
+ 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 */
}
@@ -987,16 +1010,13 @@
E->TMass[m][node] = 0.0;
for(e=1;e<=E->lmesh.nel;e++) {
- get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,
- sphere_key,rtf,E->mesh.levmax,m);
-
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]
- * dOmega.vpt[nint]
+ * E->gDA[m][e].vpt[nint]
* g_point[nint].weight[E->mesh.nsd-1]
* E->N.vpt[GNVINDEX(node,nint)];
}
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -178,7 +178,7 @@
float** SXY, float** SXZ, float** SZY,
float** divv, float** vorv)
{
- void get_global_shape_fn();
+ void get_rtf_vpts();
void velo_from_element();
void stress_conform_bcs();
@@ -189,9 +189,8 @@
double pre[9],tww[9],rtf[4][9];
double velo_scaling, stress_scaling;
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
+ struct Shape_function_dA *dOmega;
+ struct Shape_function_dx *GNx;
const int dims=E->mesh.nsd;
const int vpts=vpoints[dims];
@@ -209,9 +208,11 @@
Szy = 0.0;
div = 0.0;
vor = 0.0;
- get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,sphere_key,rtf,E->mesh.levmax,m);
+ get_rtf_vpts(E, m, lev, e, rtf);
velo_from_element(E,VV,m,e,sphere_key);
+ dOmega = &(E->gDA[m][e]);
+ GNx = &(E->gNX[m][e]);
/* Vxyz is the strain rate vector, whose relationship with
* the strain rate tensor (e) is that:
@@ -224,7 +225,7 @@
* where 1 is theta, 2 is phi, and 3 is r
*/
for(j=1;j<=vpts;j++) {
- pre[j] = E->EVi[m][(e-1)*vpts+j]*dOmega.vpt[j];
+ pre[j] = E->EVi[m][(e-1)*vpts+j]*dOmega->vpt[j];
dilation[j] = 0.0;
Vxyz[1][j] = 0.0;
Vxyz[2][j] = 0.0;
@@ -239,34 +240,34 @@
for(i=1;i<=ends;i++) {
tww[i] = 0.0;
for(j=1;j<=vpts;j++)
- tww[i] += dOmega.vpt[j] * g_point[j].weight[E->mesh.nsd-1]
+ tww[i] += dOmega->vpt[j] * g_point[j].weight[E->mesh.nsd-1]
* E->N.vpt[GNVINDEX(i,j)];
}
for(j=1;j<=vpts;j++) {
for(i=1;i<=ends;i++) {
- Vxyz[1][j]+=( VV[1][i]*GNx.vpt[GNVXINDEX(0,i,j)]
+ Vxyz[1][j]+=( VV[1][i]*GNx->vpt[GNVXINDEX(0,i,j)]
+ VV[3][i]*E->N.vpt[GNVINDEX(i,j)] )*rtf[3][j];
- Vxyz[2][j]+=( (VV[2][i]*GNx.vpt[GNVXINDEX(1,i,j)]
+ Vxyz[2][j]+=( (VV[2][i]*GNx->vpt[GNVXINDEX(1,i,j)]
+ VV[1][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j]))/sin(rtf[1][j])
+ VV[3][i]*E->N.vpt[GNVINDEX(i,j)] )*rtf[3][j];
- Vxyz[3][j]+= VV[3][i]*GNx.vpt[GNVXINDEX(2,i,j)];
+ Vxyz[3][j]+= VV[3][i]*GNx->vpt[GNVXINDEX(2,i,j)];
- Vxyz[4][j]+=( (VV[1][i]*GNx.vpt[GNVXINDEX(1,i,j)]
+ Vxyz[4][j]+=( (VV[1][i]*GNx->vpt[GNVXINDEX(1,i,j)]
- VV[2][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j]))/sin(rtf[1][j])
- + VV[2][i]*GNx.vpt[GNVXINDEX(0,i,j)])*rtf[3][j];
- Vxyz[5][j]+=VV[1][i]*GNx.vpt[GNVXINDEX(2,i,j)] + rtf[3][j]*(VV[3][i]
- *GNx.vpt[GNVXINDEX(0,i,j)]-VV[1][i]*E->N.vpt[GNVINDEX(i,j)]);
- Vxyz[6][j]+=VV[2][i]*GNx.vpt[GNVXINDEX(2,i,j)] + rtf[3][j]*(VV[3][i]
- *GNx.vpt[GNVXINDEX(1,i,j)]/sin(rtf[1][j])-VV[2][i]*E->N.vpt[GNVINDEX(i,j)]);
+ + VV[2][i]*GNx->vpt[GNVXINDEX(0,i,j)])*rtf[3][j];
+ Vxyz[5][j]+=VV[1][i]*GNx->vpt[GNVXINDEX(2,i,j)] + rtf[3][j]*(VV[3][i]
+ *GNx->vpt[GNVXINDEX(0,i,j)]-VV[1][i]*E->N.vpt[GNVINDEX(i,j)]);
+ Vxyz[6][j]+=VV[2][i]*GNx->vpt[GNVXINDEX(2,i,j)] + rtf[3][j]*(VV[3][i]
+ *GNx->vpt[GNVXINDEX(1,i,j)]/sin(rtf[1][j])-VV[2][i]*E->N.vpt[GNVINDEX(i,j)]);
Vxyz[7][j]+=rtf[3][j] * (
- VV[1][i]*GNx.vpt[GNVXINDEX(0,i,j)]
+ VV[1][i]*GNx->vpt[GNVXINDEX(0,i,j)]
+ VV[1][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j])/sin(rtf[1][j])
- + VV[2][i]*GNx.vpt[GNVXINDEX(1,i,j)]/sin(rtf[1][j]) );
+ + VV[2][i]*GNx->vpt[GNVXINDEX(1,i,j)]/sin(rtf[1][j]) );
Vxyz[8][j]+=rtf[3][j]/sin(rtf[1][j])*
- ( VV[2][i]*GNx.vpt[GNVXINDEX(0,i,j)]*sin(rtf[1][j])
+ ( VV[2][i]*GNx->vpt[GNVXINDEX(0,i,j)]*sin(rtf[1][j])
+ VV[2][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j])
- - VV[1][i]*GNx.vpt[GNVXINDEX(1,i,j)] );
+ - VV[1][i]*GNx->vpt[GNVXINDEX(1,i,j)] );
}
}
@@ -282,8 +283,8 @@
Sxy += pre[j] * Vxyz[4][j];
Sxz += pre[j] * Vxyz[5][j];
Szy += pre[j] * Vxyz[6][j];
- div += Vxyz[7][j]*dOmega.vpt[j];
- vor += Vxyz[8][j]*dOmega.vpt[j];
+ div += Vxyz[7][j]*dOmega->vpt[j];
+ vor += Vxyz[8][j]*dOmega->vpt[j];
}
Sxx /= E->eco[m][e].area;
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2008-01-22 19:38:43 UTC (rev 9110)
@@ -719,14 +719,12 @@
float *EEDOT;
int m,SQRT;
{
- void get_global_shape_fn();
+ void get_rtf_ppts();
void velo_from_element();
void construct_c3x3matrix_el();
void get_ba_p();
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
+ struct Shape_function_dx *GNx;
double edot[4][4], rtf[4][9];
double theta;
@@ -745,11 +743,9 @@
for(e=1; e<=nel; e++) {
- /* get shape function on presure gauss points */
- get_global_shape_fn(E, e, &GN, &GNx, &dOmega, 2,
- sphere_key, rtf, lev, m);
-
+ get_rtf_ppts(E, m, lev, e, rtf);
velo_from_element(E, VV, m, e, sphere_key);
+ GNx = &(E->gNX[m][e]);
theta = rtf[1][1];
@@ -782,7 +778,7 @@
construct_c3x3matrix_el(E,e,&E->element_Cc,&E->element_Ccx,lev,m,1);
}
- get_ba_p(&(E->N), &GNx, &E->element_Cc, &E->element_Ccx,
+ get_ba_p(&(E->N), GNx, &E->element_Cc, &E->element_Ccx,
rtf, E->mesh.nsd, ba);
for(j=1;j<=ppts;j++)
@@ -796,27 +792,27 @@
else {
for(j=1; j<=ppts; j++) {
for(i=1; i<=ends; i++) {
- Vxyz[1][j] += (VV[1][i] * GNx.ppt[GNPXINDEX(0, i, j)]
+ Vxyz[1][j] += (VV[1][i] * GNx->ppt[GNPXINDEX(0, i, j)]
+ VV[3][i] * E->N.ppt[GNPINDEX(i, j)])
* rtf[3][j];
- Vxyz[2][j] += ((VV[2][i] * GNx.ppt[GNPXINDEX(1, i, j)]
+ Vxyz[2][j] += ((VV[2][i] * GNx->ppt[GNPXINDEX(1, i, j)]
+ VV[1][i] * E->N.ppt[GNPINDEX(i, j)]
* cos(rtf[1][j])) / sin(rtf[1][j])
+ VV[3][i] * E->N.ppt[GNPINDEX(i, j)])
* rtf[3][j];
- Vxyz[3][j] += VV[3][i] * GNx.ppt[GNPXINDEX(2, i, j)];
+ Vxyz[3][j] += VV[3][i] * GNx->ppt[GNPXINDEX(2, i, j)];
- Vxyz[4][j] += ((VV[1][i] * GNx.ppt[GNPXINDEX(1, i, j)]
+ Vxyz[4][j] += ((VV[1][i] * GNx->ppt[GNPXINDEX(1, i, j)]
- VV[2][i] * E->N.ppt[GNPINDEX(i, j)]
* cos(rtf[1][j])) / sin(rtf[1][j])
- + VV[2][i] * GNx.ppt[GNPXINDEX(0, i, j)])
+ + VV[2][i] * GNx->ppt[GNPXINDEX(0, i, j)])
* rtf[3][j];
- Vxyz[5][j] += VV[1][i] * GNx.ppt[GNPXINDEX(2, i, j)]
- + rtf[3][j] * (VV[3][i] * GNx.ppt[GNPXINDEX(0, i, j)]
+ Vxyz[5][j] += VV[1][i] * GNx->ppt[GNPXINDEX(2, i, j)]
+ + rtf[3][j] * (VV[3][i] * GNx->ppt[GNPXINDEX(0, i, j)]
- VV[1][i] * E->N.ppt[GNPINDEX(i, j)]);
- Vxyz[6][j] += VV[2][i] * GNx.ppt[GNPXINDEX(2, i, j)]
+ Vxyz[6][j] += VV[2][i] * GNx->ppt[GNPXINDEX(2, i, j)]
+ rtf[3][j] * (VV[3][i]
- * GNx.ppt[GNPXINDEX(1, i, j)]
+ * GNx->ppt[GNPXINDEX(1, i, j)]
/ sin(rtf[1][j])
- VV[2][i] * E->N.ppt[GNPINDEX(i, j)]);
}
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2008-01-22 17:06:50 UTC (rev 9109)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2008-01-22 19:38:43 UTC (rev 9110)
@@ -768,6 +768,11 @@
float *age[NCS]; /* nodal weightings */
float *age_t;
+ struct Shape_function_dx *GNX[MAX_LEVELS][NCS];
+ struct Shape_function_dA *GDA[MAX_LEVELS][NCS];
+ struct Shape_function_dx *gNX[NCS];
+ struct Shape_function_dA *gDA[NCS];
+
struct Shape_function1 M; /* master-element shape funtions */
struct Shape_function1_dx Mx;
struct Shape_function N;
More information about the cig-commits
mailing list