[cig-commits] r15473 - mc/3D/CitcomS/trunk/lib

[becker at geodynamics.org]

Author: becker
Date: 2009-07-23 11:40:41 -0700 (Thu, 23 Jul 2009)
New Revision: 15473

Modified:
   mc/3D/CitcomS/trunk/lib/Topo_gravity.c
   mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
Log:
Corrected call to get_ba in Topo_gravity for precise_strain_rate, still testing.
Rest should be cosmetic.



Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c	2009-07-23 03:50:26 UTC (rev 15472)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c	2009-07-23 18:40:41 UTC (rev 15473)
@@ -215,6 +215,8 @@
   void get_rtf_at_vpts();
   void velo_from_element();
   void stress_conform_bcs();
+  void construct_c3x3matrix_el();
+  void get_ba();
 
   int i,j,p,q,e,node,m;
 
@@ -228,15 +230,17 @@
   struct Shape_function_dA *dOmega;
   struct Shape_function_dx *GNx;
 
-  const int dims=E->mesh.nsd;
-  const int vpts=vpoints[dims];
-  const int ends=enodes[dims];
-  const int lev=E->mesh.levmax;
+
+
   const int sphere_key=1;
-  const int ppts = ppoints[dims];
+  const int lev = E->mesh.levmax;
+  const int dims = E->mesh.nsd;
+  const int nel = E->lmesh.nel;
+  const int vpts = vpoints[dims];
+  const int ends = enodes[dims];
 
   for(m=1;m<=E->sphere.caps_per_proc;m++) {
-    for(e=1;e<=E->lmesh.nel;e++)  {
+    for(e=1;e <= nel;e++)  {
       Szz = 0.0;
       Sxx = 0.0;
       Syy = 0.0;
@@ -265,7 +269,8 @@
        *    Vxyz[6] = 2*e23
        * where 1 is theta, 2 is phi, and 3 is r
        */
-      for(j=1;j<=vpts;j++)  {
+      for(j=1;j <= vpts;j++)  {	/* loop through velocity Gauss points  */
+
         pre[j] =  E->EVi[m][(e-1)*vpts+j]*dOmega->vpt[j];
         dilation[j] = 0.0;
         Vxyz[1][j] = 0.0;
@@ -280,7 +285,7 @@
 
       for(i=1;i<=ends;i++) {
         tww[i] = 0.0;
-        for(j=1;j<=vpts;j++)	/* weighting, consisting of Jacobian,
+        for(j=1;j <= vpts;j++)	/* weighting, consisting of Jacobian,
 				   Gauss weight and shape function,
 				   evaluated at integration points */
           tww[i] += dOmega->vpt[j] * g_point[j].weight[E->mesh.nsd-1]
@@ -288,28 +293,29 @@
       }
       if (E->control.precise_strain_rate){
 	/*  
-	    Ba method 
+	    B method 
 	*/
 
 	if ((e-1)%E->lmesh.elz==0) {
 	  construct_c3x3matrix_el(E,e,&E->element_Cc,&E->element_Ccx,lev,m,1);
 	}
+	/* get b at velocity Gauss points */
+	get_ba(&(E->N), GNx, &E->element_Cc, &E->element_Ccx,rtf, dims, ba);
 	
-	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<=8;p++)
-	    for(i=1;i<=ends;i++)
-	      for(q=1;q<=dims;q++) {
+	for(j=1;j <= vpts;j++)
+	  for(p=1;p <= 8;p++)
+	    for(i=1;i <= ends;i++)
+	      for(q=1;q <= dims;q++) {
 		Vxyz[p][j] += ba[i][j][q][p] * VV[q][i];
 	      }
 	
 	
-      }else{		/* old method */
+      }else{	
+	/* old method */
+
 	/* integrate over element  */
-	for(j=1;j<=vpts;j++)   {	/* Gauss integration points */
-	  for(i=1;i<=ends;i++)   { /* nodes in element loop */
+	for(j=1;j <= vpts;j++)   {	/* Gauss integration points */
+	  for(i=1;i <= ends;i++)   { /* nodes in element loop */
 	    /* strain rate contributions from each node */
 	    Vxyz[1][j]+=( VV[1][i]*GNx->vpt[GNVXINDEX(0,i,j)]
 			  + VV[3][i]*E->N.vpt[GNVINDEX(i,j)] )*rtf[3][j];
@@ -342,7 +348,8 @@
               dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
       }
 
-      for(j=1;j<=vpts;j++)   {
+      for(j=1;j <= vpts;j++)   {
+	/* deviatoric stress, pressure will be added later */
           Sxx += 2.0 * pre[j] * (Vxyz[1][j] - dilation[j]); /*  */
           Syy += 2.0 * pre[j] * (Vxyz[2][j] - dilation[j]);
           Szz += 2.0 * pre[j] * (Vxyz[3][j] - dilation[j]);

Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2009-07-23 03:50:26 UTC (rev 15472)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2009-07-23 18:40:41 UTC (rev 15473)
@@ -941,7 +941,8 @@
 
     for(e=1; e<=nel; e++) {
 
-        get_rtf_at_ppts(E, m, lev, e, rtf);
+        get_rtf_at_ppts(E, m, lev, e, rtf); /* pressure evaluation
+					       points */
         velo_from_element(E, VV, m, e, sphere_key);
         GNx = &(E->gNX[m][e]);
 
@@ -1019,7 +1020,7 @@
                                        - VV[2][i] * E->N.ppt[GNPINDEX(i, j)]);
                 }
             }
-        } /* end of else */
+        } /* end of fast, imprecise strain-rate computation */
 
         if(E->control.inv_gruneisen != 0) {
             for(j=1; j<=ppts; j++)
@@ -1030,7 +1031,7 @@
         edot[1][2] = edot[1][3] = edot[2][3] = 0;
 
         /* edot is 2 * (the deviatoric strain rate tensor) */
-        for(j=1; j<=ppts; j++) {
+        for(j=1; j <= ppts; j++) {
             edot[1][1] += 2.0 * (Vxyz[1][j] - dilation[j]);
             edot[2][2] += 2.0 * (Vxyz[2][j] - dilation[j]);
             edot[3][3] += 2.0 * (Vxyz[3][j] - dilation[j]);