[cig-commits] r7836 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Aug 16 13:10:58 PDT 2007
Author: tan2
Date: 2007-08-16 13:10:58 -0700 (Thu, 16 Aug 2007)
New Revision: 7836
Modified:
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
Log:
When the element is close to the poles, use a more precise method to compute the strain rate.
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-08-16 19:38:56 UTC (rev 7835)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-08-16 20:10:58 UTC (rev 7836)
@@ -224,6 +224,75 @@
/*==============================================================
+ Function to supply the element strain-displacement matrix Ba at pressure
+ quadrature points, which is used to compute strain rate
+ ============================================================== */
+
+void get_ba_p(struct Shape_function *N, struct Shape_function_dx *GNx,
+ struct CC *cc, struct CCX *ccx, double rtf[4][9],
+ int dims, double ba[9][9][4][7])
+{
+ int k, a, n;
+ const int ppts = PPOINTS3D;
+ const int ends = ENODES3D;
+
+ double ra[9], isi[9], ct[9];
+ double gnx0, gnx1, gnx2, shp, cc1, cc2, cc3;
+
+ for(k=1;k<=ppts;k++) {
+ ra[k] = rtf[3][k];
+ isi[k] = 1.0 / sin(rtf[1][k]);
+ ct[k] = cos(rtf[1][k]) * isi[k];
+ }
+
+ for(k=1;k<=ppts;k++)
+ for(a=1;a<=ends;a++) {
+ gnx0 = GNx->ppt[GNPXINDEX(0,a,k)];
+ gnx1 = GNx->ppt[GNPXINDEX(1,a,k)];
+ gnx2 = GNx->ppt[GNPXINDEX(2,a,k)];
+ shp = N->ppt[GNPINDEX(a,k)];
+ for(n=1;n<=dims;n++) {
+ cc1 = cc->ppt[BPINDEX(1,n,a,k)];
+ cc2 = cc->ppt[BPINDEX(2,n,a,k)];
+ cc3 = cc->ppt[BPINDEX(3,n,a,k)];
+
+ ba[a][k][n][1] = ( gnx0 * cc1
+ + shp * ccx->ppt[BPXINDEX(1,n,1,a,k)]
+ + shp * cc3 ) * ra[k];
+
+ ba[a][k][n][2] = ( shp * cc1 * ct[k]
+ + shp * cc3
+ + ( gnx1 * cc2
+ + shp * ccx->ppt[BPXINDEX(2,n,2,a,k)] )
+ * isi[k] ) * ra[k];
+
+ ba[a][k][n][3] = gnx2 * cc3;
+
+ ba[a][k][n][4] = ( gnx0 * cc2
+ + shp * ccx->ppt[BPXINDEX(2,n,1,a,k)]
+ - shp * cc2 * ct[k]
+ + ( gnx1 * cc1
+ + shp * ccx->ppt[BPXINDEX(1,n,2,a,k)] )
+ * isi[k] ) * ra[k];
+
+ ba[a][k][n][5] = gnx2 * cc1
+ + ( gnx0 * cc3
+ + shp * ( ccx->ppt[BPXINDEX(3,n,1,a,k)]
+ - cc1 ) ) * ra[k];
+
+ ba[a][k][n][6] = gnx2 * cc2
+ - ra[k] * shp * cc2
+ + ( gnx1 * cc3
+ + shp * ccx->ppt[BPXINDEX(3,n,2,a,k)] )
+ * isi[k] * ra[k];
+ }
+ }
+ return;
+}
+
+
+
+/*==============================================================
Function to supply the element k matrix for a given element e.
============================================================== */
@@ -249,8 +318,6 @@
struct Shape_function_dx GNx;
double ba[9][9][4][7]; /* integration points,node,3x6 matrix */
- /* d1 and d2 are the elastic coefficient matrix for incompressible
- and compressible creeping flow, respectively */
const int nn=loc_mat_size[E->mesh.nsd];
const int vpts = VPOINTS3D;
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-16 19:38:56 UTC (rev 7835)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-16 20:10:58 UTC (rev 7836)
@@ -702,14 +702,15 @@
{
void get_global_shape_fn();
void velo_from_element();
+ void get_ba_p();
struct Shape_function GN;
struct Shape_function_dA dOmega;
struct Shape_function_dx GNx;
-
-
double edot[4][4], rtf[4][9];
+ double theta;
+ double ba[9][9][4][7];
float VV[4][9], Vxyz[7][9], dilation[9];
int e, i, j, p, q, n;
@@ -730,7 +731,9 @@
velo_from_element(E, VV, m, e, sphere_key);
+ theta = rtf[1][1];
+
/* Vxyz is the strain rate vector, whose relationship with
* the strain rate tensor (e) is that:
* Vxyz[1] = e11
@@ -751,34 +754,55 @@
dilation[j] = 0.0;
}
- for(j=1; j<=ppts; j++) {
- for(i=1; i<=ends; i++) {
- 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)]
- + 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)];
+ if ((theta < 0.09) || (theta > 3.05)) {
+ /* When the element is close to the poles, use a more
+ * precise method to compute the strain rate. */
- 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)])
- * 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)]
- - VV[1][i] * E->N.ppt[GNPINDEX(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)]
- / sin(rtf[1][j])
- - VV[2][i] * E->N.ppt[GNPINDEX(i, j)]);
+ if ((e-1)%E->lmesh.ELZ[lev]==0) {
+ 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,
+ rtf, E->mesh.nsd, ba);
+
+ for(j=1;j<=ppts;j++)
+ for(p=1;p<=6;p++)
+ for(i=1;i<=ends;i++)
+ for(q=1;q<=dims;q++) {
+ Vxyz[p][j] += ba[i][j][p][q] * VV[q][i];
+ }
+
}
+ 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)]
+ + VV[3][i] * E->N.ppt[GNPINDEX(i, j)])
+ * rtf[3][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[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)])
+ * 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)]
+ - VV[1][i] * E->N.ppt[GNPINDEX(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)]
+ / sin(rtf[1][j])
+ - VV[2][i] * E->N.ppt[GNPINDEX(i, j)]);
+ }
+ }
+ } /* end of else */
+
if(E->control.inv_gruneisen != 0) {
for(j=1; j<=ppts; j++)
dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
@@ -797,7 +821,6 @@
edot[2][3] += Vxyz[6][j];
}
- /* TODO: figure out why 2.0 factor is here */
EEDOT[e] = edot[1][1] * edot[1][1]
+ edot[1][2] * edot[1][2] * 2.0
+ edot[2][2] * edot[2][2]
More information about the cig-commits
mailing list