/*
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 *
 *<LicenseText>
 *
 * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
 * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
 * Copyright (C) 1994-2005, California Institute of Technology.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *</LicenseText>
 *
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 */
/* Functions relating to the determination of viscosity field either
   as a function of the run, as an initial condition or as specified from
   a previous file */


#include <math.h>
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
#include "parsing.h"


void myerror(struct All_variables *,char *);

void allocate_visc_vars(struct All_variables *E);
void read_visc_layer_file(struct All_variables *E);
void read_visc_param_from_file(struct All_variables *E,
                               const char *param, float *var,
                               FILE *fp);
float compute_nodal_stress_state_factor(struct All_variables *,int,int);
static void apply_low_visc_wedge_channel(struct All_variables *E, float **evisc);
static void low_viscosity_channel_factor(struct All_variables *E, float *F);
static void low_viscosity_wedge_factor(struct All_variables *E, float *F);
void parallel_process_termination();
#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
#include "anisotropic_viscosity.h"
#endif


void viscosity_system_input(struct All_variables *E)
{
    int m=E->parallel.me;
    int i;


    input_boolean("VISC_UPDATE",&(E->viscosity.update_allowed),"on",m);

    input_boolean("visc_layer_control",&(E->viscosity.layer_control),"off",m);
    input_string("visc_layer_file", E->viscosity.layer_file,"visc.dat",m);


    input_int("allow_anisotropic_viscosity",&(E->viscosity.allow_anisotropic_viscosity),"0",m);
#ifndef CITCOM_ALLOW_ANISOTROPIC_VISC 
    if(E->viscosity.allow_anisotropic_viscosity){ /* error */
      fprintf(stderr,"error: allow_anisotropic_viscosity is not zero, but code not compiled with CITCOM_ALLOW_ANISOTROPIC_VISC\n");
      parallel_process_termination();
    }
#else
    if(E->viscosity.allow_anisotropic_viscosity){ /* read additional
						     parameters for
						     anisotropic
						     viscosity */
      input_int("anisotropic_init",&(E->viscosity.anisotropic_init),"0",m); /* 0: isotropic
									       1: random
									       2: read in director orientation
									       and log10(eta_s/eta) 
									       3: align with velocity, use ani_vis2_factor for eta_s/eta
									       4: align with ISA, use ani_vis2_factor for eta_s/eta
									       5: align mixed depending on deformation state, use ani_vis2_factor for eta_s/eta
									       6: use radially aligned director and taper eta_s/eta from base (1) to top of layer (ani_vis2_factor)
									       
									    */
      input_string("anisotropic_init_dir",(E->viscosity.anisotropic_init_dir),"",m); /* directory
											for
											ggrd
											type
											init */
      input_int("anivisc_layer",&(E->viscosity.anivisc_layer),"1",m); /* >0: assign to layers on top of anivisc_layer
									 <0: assign to layer = anivisc_layer
								      */
      if((E->viscosity.anisotropic_init == 6) && (E->viscosity.anivisc_layer >= 0))
	myerror(E,"anisotropic init mode 6 requires selection of layer where anisotropy applies");
      
      input_boolean("anivisc_start_from_iso",
		    &(E->viscosity.anivisc_start_from_iso),"on",m); /* start
								       from
								       isotropic
								       solution? */
      if(!E->viscosity.anivisc_start_from_iso)
	if(E->viscosity.anisotropic_init == 3){
	  if(E->parallel.me == 0)fprintf(stderr,"WARNING: overriding anisotropic first step for ani init mode\n");
	  E->viscosity.anivisc_start_from_iso = TRUE;
	}
      /* ratio between weak and strong direction */
      input_double("ani_vis2_factor",&(E->viscosity.ani_vis2_factor),"1.0",m);
      
      
    }
#endif
    /* allocate memory here */
    allocate_visc_vars(E);

    for(i=0;i < E->viscosity.num_mat;i++)
        E->viscosity.N0[i]=1.0;
    input_float_vector("visc0",E->viscosity.num_mat,(E->viscosity.N0),m);

    input_boolean("TDEPV",&(E->viscosity.TDEPV),"on",m);
    input_int("rheol",&(E->viscosity.RHEOL),"3",m);
    if (E->viscosity.TDEPV) {
        for(i=0;i < E->viscosity.num_mat;i++) {
            E->viscosity.T[i] = 0.0;
            E->viscosity.Z[i] = 0.0;
            E->viscosity.E[i] = 0.0;
	    E->viscosity.Te[i] = 0.0;
        }

        input_float_vector("viscT",E->viscosity.num_mat,(E->viscosity.T),m);
        input_float_vector("viscE",E->viscosity.num_mat,(E->viscosity.E),m);
        input_float_vector("viscZ",E->viscosity.num_mat,(E->viscosity.Z),m);
        input_float_vector("viscTe",E->viscosity.num_mat,(E->viscosity.Te),m);

        /* for viscosity 8 */
        input_float("T_sol0",&(E->viscosity.T_sol0),"0.6",m);
        input_float("ET_red",&(E->viscosity.ET_red),"0.1",m);
    }


    E->viscosity.sdepv_misfit = 1.0;
    input_boolean("SDEPV",&(E->viscosity.SDEPV),"off",m);
    if (E->viscosity.SDEPV) {
      E->viscosity.sdepv_visited = 0;
      input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
      input_int("sdepv_rheol",&(E->viscosity.sdepv_rheol),"1",m); /* different types of stress dependence 
								     1: default (for backward compatibility/doesn't make sense?) 
								     2: new implementation */
      input_float_vector("sdepv_trns",E->viscosity.num_mat,(E->viscosity.sdepv_trns),m); /* transition strain-rate */
      input_boolean("sdepv_start_from_newtonian",&(E->viscosity.sdepv_start_from_newtonian),"on",m);
    }
    input_boolean("visc_adjust_N0",&(E->viscosity.adjust_N0),"off",m); /* iterate on prefactors */
    input_float("adjust_N0_alpha",&(E->viscosity.adjust_N0_alpha),"1.0",m); /* damping factor, not sure 
									       if this should be simply unity, or maybe something like 0.95 */
    
    input_boolean("PDEPV",&(E->viscosity.PDEPV),"off",m); /* plasticity addition by TWB */
    /* Choose weakening, i.e., 0 - off; 1 - PSS; 2 - VSS; 3 - VSS^power */
    input_int("strain_dep_plasticity",&(E->viscosity.strain_dep_plasticity),"0",m);
    if (E->viscosity.PDEPV) {
      /* strain dependent plasticity addition by LF */
      /* The following two lines have been changed!*/
      input_float("gamma_cr",&(E->viscosity.gamma_cr),"10",m);
      input_float("maxD",&(E->viscosity.maxD),"0.9",m); 
      E->viscosity.pdepv_visited = 0;
      for(i=0;i < E->viscosity.num_mat;i++) {
          E->viscosity.pdepv_a[i] = 1.e20; /* \sigma_y = min(a + b * (1-r),y) */
          E->viscosity.pdepv_b[i] = 0.0;
          E->viscosity.pdepv_y[i] = 1.e20;
      }

      input_boolean("psrw",&(E->viscosity.psrw),"off",m); /* SRW? else regular plasiticity */
      input_boolean("pdepv_eff",&(E->viscosity.pdepv_eff),"on",m); /* effective
								      or
								      min/max
								      approach */
      input_float_vector("pdepv_a",E->viscosity.num_mat,(E->viscosity.pdepv_a),m);
      input_float_vector("pdepv_b",E->viscosity.num_mat,(E->viscosity.pdepv_b),m);
      input_float_vector("pdepv_y",E->viscosity.num_mat,(E->viscosity.pdepv_y),m);

      input_float("pdepv_offset",&(E->viscosity.pdepv_offset),"0.0",m);

      input_boolean("higher_cont_yield",&(E->control.stiff_cont),"off",m);
      input_float("cont_sig_y",&(E->viscosity.sig_c),"1e7",m);
    }

    /* stress state dependent weakening? 1: for damage 2: for plasticity */
    input_int("p_stress_state_dependent",&(E->viscosity.p_stress_state_dependent),"0",m); 

    input_float("sdepv_misfit",&(E->viscosity.sdepv_misfit),"0.001",m);	/* there should be no harm in having 
									   this parameter read in regardless of 
									   rheology (activated it for anisotropic viscosity)
									*/

    // moved to Composition related, for init purposes
    input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
    for(i=0;i<10;i++)
      E->viscosity.cdepv_ff[i] = 1.0; /* flavor factors for CDEPV */
    if(E->viscosity.CDEPV){
      /* compositional viscosity */
      if(E->control.tracer < 1){
	fprintf(stderr,"error: CDEPV requires tracers, but tracer is off\n");
	parallel_process_termination();
      }
      if(E->trace.nflavors > 10)
	myerror(E,"error: too many flavors for CDEPV");
      /* read in flavor factors */
      input_float_vector("cdepv_ff",E->trace.nflavors,
			 (E->viscosity.cdepv_ff),m);
      /* and take the log because we're using a geometric avg */
      for(i=0;i<E->trace.nflavors;i++)
	E->viscosity.cdepv_ff[i] = log(E->viscosity.cdepv_ff[i]);
    }


    input_boolean("low_visc_channel",&(E->viscosity.channel),"off",m);
    input_boolean("low_visc_wedge",&(E->viscosity.wedge),"off",m);

    input_float("lv_min_radius",&(E->viscosity.lv_min_radius),"0.9764",m);
    input_float("lv_max_radius",&(E->viscosity.lv_max_radius),"0.9921",m);
    input_float("lv_channel_thickness",&(E->viscosity.lv_channel_thickness),"0.0047",m);
    input_float("lv_reduction",&(E->viscosity.lv_reduction),"0.5",m);

    input_boolean("use_ne_visc_smooth",&(E->viscosity.SMOOTH),"off",m);

    input_boolean("VMAX",&(E->viscosity.MAX),"off",m);
    if (E->viscosity.MAX)
        input_float("visc_max",&(E->viscosity.max_value),"1e22,1,nomax",m);

    input_boolean("VMIN",&(E->viscosity.MIN),"off",m);
    if (E->viscosity.MIN)
        input_float("visc_min",&(E->viscosity.min_value),"1e20",m);



    
    
    return;
}


void viscosity_input(struct All_variables *E)
{
    int m = E->parallel.me;

    input_string("Viscosity",E->viscosity.STRUCTURE,"system",m);
    input_int ("visc_smooth_method",&(E->viscosity.smooth_cycles),"0",m);

    if ( strcmp(E->viscosity.STRUCTURE,"system") == 0)
        E->viscosity.FROM_SYSTEM = 1;
    else
        E->viscosity.FROM_SYSTEM = 0;

    if (E->viscosity.FROM_SYSTEM)
        viscosity_system_input(E);

    return;
}


void allocate_visc_vars(struct All_variables *E)
{
  int i,j,k,lim,nel,nno;

  if(E->viscosity.layer_control) {
    /* noz is already defined, but elz is undefined yet */
    i = E->mesh.noz - 1;
  } else {
    i = E->viscosity.num_mat;
    
  }
  i++;
  
  E->viscosity.N0 = (float*) malloc(i*sizeof(float));
  E->viscosity.E = (float*) malloc(i*sizeof(float));
  E->viscosity.T = (float*) malloc(i*sizeof(float));
  E->viscosity.Te = (float*) malloc(i*sizeof(float));
  E->viscosity.Z = (float*) malloc(i*sizeof(float));
  E->viscosity.pdepv_a = (float*) malloc(i*sizeof(float));
  E->viscosity.pdepv_b = (float*) malloc(i*sizeof(float));
  E->viscosity.pdepv_y = (float*) malloc(i*sizeof(float));
  E->viscosity.sdepv_expt = (float*) malloc(i*sizeof(float));
  E->viscosity.sdepv_trns = (float*) malloc(i*sizeof(float));
  
  if(E->viscosity.N0 == NULL ||
     E->viscosity.E == NULL ||
     E->viscosity.T == NULL ||
     E->viscosity.Te == NULL ||
     E->viscosity.Z == NULL ||
     E->viscosity.pdepv_a == NULL ||
     E->viscosity.pdepv_b == NULL ||
     E->viscosity.pdepv_y == NULL ||
     E->viscosity.sdepv_expt == NULL) {
    fprintf(stderr, "Error: Cannot allocate visc memory, rank=%d\n",
	    E->parallel.me);
    parallel_process_termination();
  }

}


/* ============================================ */
/* 

   when called from Drive_solvers:
   
   evisc = E->EVI[E->mesh.levmax]
   visc  = E->VI[E->mesh.levmax]

 */
void get_system_viscosity(E,propogate,evisc,visc)
     struct All_variables *E;
     int propogate;
     float **evisc,**visc;
{
    void visc_from_mat();
    void visc_from_T();
    void visc_from_S();

    void visc_from_P();
    void visc_from_C();

    void apply_viscosity_smoother();
    void visc_from_gint_to_nodes();
    void visc_from_nodes_to_gint();
    void viscy_from_element_to_nodes();

    int i,j,m;
    float temp1,temp2,*vvvis;
    double *TG;
    float *SXX[NCS],*SYY[NCS],*SXY[NCS],*SXZ[NCS],*SZY[NCS],*SZZ[NCS];
    float *divv[NCS],*vorv[NCS];

    const int vpts = vpoints[E->mesh.nsd];
#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
    if(E->viscosity.allow_anisotropic_viscosity){
      if(!E->viscosity.anisotropic_viscosity_init)
	set_anisotropic_viscosity_at_element_level(E,1);
      else
	set_anisotropic_viscosity_at_element_level(E,0);
    }
#endif


    if(E->viscosity.TDEPV)
        visc_from_T(E,evisc,propogate);
    else
        visc_from_mat(E,evisc);

    if(E->viscosity.CDEPV)	/* compositional prefactor */
      visc_from_C(E,evisc);

    if(E->viscosity.SDEPV)
      visc_from_S(E,evisc,propogate);

    if(E->viscosity.PDEPV)	/* "plasticity" */
      visc_from_P(E,evisc);


    /* i think this should me placed differently i.e.  before the
       stress dependence but I won't change it because it's by
       someone else

       TWB
    */
    if(E->viscosity.channel || E->viscosity.wedge)
        apply_low_visc_wedge_channel(E, evisc);


    /* min/max cut-off */

    if(E->viscosity.MAX) {
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=E->lmesh.nel;i++)
                for(j=1;j<=vpts;j++)
                    if(evisc[m][(i-1)*vpts + j] > E->viscosity.max_value)
                        evisc[m][(i-1)*vpts + j] = E->viscosity.max_value;
    }

    if(E->viscosity.MIN) {
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=E->lmesh.nel;i++)
                for(j=1;j<=vpts;j++)
                    if(evisc[m][(i-1)*vpts + j] < E->viscosity.min_value)
                        evisc[m][(i-1)*vpts + j] = E->viscosity.min_value;
    }
    if(E->viscosity.adjust_N0 &&
       (E->monitor.solution_cycles == E->monitor.solution_cycles_init))
      adjust_visc_N0(E,evisc);	/* adjust N0 to match layers on next
				   step, for first time step */

    if (E->control.verbose)  {
      fprintf(E->fp_out,"output_evisc \n");
      for(m=1;m<=E->sphere.caps_per_proc;m++) {
        fprintf(E->fp_out,"output_evisc for cap %d\n",E->sphere.capid[m]);
      for(i=1;i<=E->lmesh.nel;i++)
          fprintf(E->fp_out,"%d %d %f %f\n",i,E->mat[m][i],evisc[m][(i-1)*vpts+1],evisc[m][(i-1)*vpts+7]);
      }
      fflush(E->fp_out);
    }
    /* interpolate from gauss quadrature points to node points for output */
    visc_from_gint_to_nodes(E,evisc,visc,E->mesh.levmax);

    if(E->viscosity.SMOOTH){ /* go the other way, for
					    smoothing */
      visc_from_nodes_to_gint(E,visc,evisc,E->mesh.levmax);
    }

#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC /* allow for anisotropy */
    if(E->viscosity.allow_anisotropic_viscosity){
      visc_from_gint_to_nodes(E,E->EVI2[E->mesh.levmax], E->VI2[E->mesh.levmax],E->mesh.levmax);
      visc_from_gint_to_nodes(E,E->EVIn1[E->mesh.levmax], E->VIn1[E->mesh.levmax],E->mesh.levmax);
      visc_from_gint_to_nodes(E,E->EVIn2[E->mesh.levmax], E->VIn2[E->mesh.levmax],E->mesh.levmax);
      visc_from_gint_to_nodes(E,E->EVIn3[E->mesh.levmax], E->VIn3[E->mesh.levmax],E->mesh.levmax);
      normalize_director_at_nodes(E,E->VIn1[E->mesh.levmax],E->VIn2[E->mesh.levmax],E->VIn3[E->mesh.levmax],E->mesh.levmax);
      
      if(E->viscosity.SMOOTH){ 
	if(E->parallel.me == 0)fprintf(stderr,"WARNING: smoothing anisotropic viscosity, perhaps not a good idea\n");
	visc_from_nodes_to_gint(E,E->VI2[E->mesh.levmax],E->EVI2[E->mesh.levmax],E->mesh.levmax);
	visc_from_nodes_to_gint(E,E->VIn1[E->mesh.levmax],E->EVIn1[E->mesh.levmax],E->mesh.levmax);
	visc_from_nodes_to_gint(E,E->VIn2[E->mesh.levmax],E->EVIn2[E->mesh.levmax],E->mesh.levmax);
	visc_from_nodes_to_gint(E,E->VIn3[E->mesh.levmax],E->EVIn3[E->mesh.levmax],E->mesh.levmax);
	normalize_director_at_gint(E,E->EVIn1[E->mesh.levmax],E->EVIn2[E->mesh.levmax],E->EVIn3[E->mesh.levmax],E->mesh.levmax);
    
      }
    }
#endif

    if(E->viscosity.p_stress_state_dependent){
      /* compute nodal stress for next time around */
      allocate_STD_mem(E, SXX, SYY, SZZ, SXY, SXZ, SZY, divv, vorv);
      compute_nodal_stress(E, SXX, SYY, SZZ, SXY, SXZ, SZY, divv, vorv);
      free_STD_mem(E, SXX, SYY, SZZ, SXY, SXZ, SZY, divv, vorv);
    }
    return;
}



void initial_viscosity(struct All_variables *E)
{
    void report(struct All_variables*, char*);

    report(E,"Initialize viscosity field");

    if (E->viscosity.FROM_SYSTEM)
        get_system_viscosity(E,1,E->EVI[E->mesh.levmax],E->VI[E->mesh.levmax]);

    return;
}


void visc_from_mat(E,EEta)
     struct All_variables *E;
     float **EEta;
{

    int i,m,jj;
    if(E->control.mat_control){	/* use pre-factor even without temperature dependent viscosity */
      for(m=1;m<=E->sphere.caps_per_proc;m++)
        for(i=1;i<=E->lmesh.nel;i++)
	  for(jj=1;jj<=vpoints[E->mesh.nsd];jj++)
	    EEta[m][ (i-1)*vpoints[E->mesh.nsd]+jj ] = E->viscosity.N0[E->mat[m][i]-1]*E->VIP[m][i];
     }else{
      for(m=1;m<=E->sphere.caps_per_proc;m++)
        for(i=1;i<=E->lmesh.nel;i++)
	  for(jj=1;jj<=vpoints[E->mesh.nsd];jj++)
	    EEta[m][ (i-1)*vpoints[E->mesh.nsd]+jj ] = E->viscosity.N0[E->mat[m][i]-1];
    }

    return;
}
/* 
   
   adjust N0 and viscosities to match original layer N0 on average
   
*/
void adjust_visc_N0(E,EEta)
     struct All_variables *E;
     float **EEta;
{
  int ie,m,jj,ilay;
  float laya[CITCOM_MAX_VISC_LAYER],glaya[CITCOM_MAX_VISC_LAYER],scale;
  int laya_c[CITCOM_MAX_VISC_LAYER],glaya_c[CITCOM_MAX_VISC_LAYER];
  static int init=0;
  static float N0_orig[CITCOM_MAX_VISC_LAYER];

  const int vp = vpoints[E->mesh.nsd];
  const int nel = E->lmesh.nel;

  /* compute average for all layers */
  for(ilay=0;ilay < E->viscosity.num_mat;ilay++){
    laya[ilay]  = glaya[ilay] = 0.;
    laya_c[ilay] = glaya_c[ilay] = 0;
    if(!init)
      N0_orig[ilay] = E->viscosity.N0[ilay]; /* store original prefactor */
  }
  for(m=1;m<=E->sphere.caps_per_proc;m++){
    for(ie=1;ie <= nel;ie++){
      ilay = E->mat[m][ie]-1;
      for(jj=1;jj <= vp;jj++){
	laya[ilay] += log(EEta[m][(ie-1)*vp+jj]); /* log average */
	laya_c[ilay]++;
      }
    }
  }
  MPI_Allreduce(laya,glaya,E->viscosity.num_mat,
		MPI_FLOAT,MPI_SUM,E->parallel.horizontal_comm);
  MPI_Allreduce(laya_c,glaya_c,E->viscosity.num_mat,
		MPI_INT,MPI_SUM,E->parallel.horizontal_comm);
  for(m=1;m<=E->sphere.caps_per_proc;m++){
    for(ilay=0;ilay < E->viscosity.num_mat;ilay++){
      if(glaya_c[ilay]){
	glaya[ilay] /= (float)glaya_c[ilay];
	glaya[ilay] = exp(glaya[ilay]);
	
	if(init > 1){/* adjust after second iteration */
	  scale = pow(10,E->viscosity.adjust_N0_alpha*log10(N0_orig[ilay]/glaya[ilay]));
	  for(ie=1;ie <= nel;ie++){	/* rescale viscosities */
	    if(E->mat[m][ie] == ilay+1){
	      for(jj=1;jj <= vp;jj++){
		EEta[m][(ie-1)*vp+jj] *= scale; /*  */
		if(E->viscosity.MIN)
		  if(EEta[m][(ie-1)*vp+jj] < E->viscosity.min_value)
		    EEta[m][(ie-1)*vp+jj] = E->viscosity.min_value;
		if(E->viscosity.MAX)
		  if(EEta[m][(ie-1)*vp+jj] > E->viscosity.max_value)
		    EEta[m][(ie-1)*vp+jj] = E->viscosity.max_value;
	      }
	    }
	  }
	  E->viscosity.N0[ilay] *= scale;
	}else{
	  scale = 1;
	}
	
	if((E->parallel.me_loc[1]==0)&&(E->parallel.me_loc[2]==0)) /* select one of the horizontal processor chunks */
	  fprintf(stderr,"adjust_visc_N0-%i: lay %i: visc_avg: %10.3e N0_orig: %8.2f N0: %8.2f scale: %.2e alpha: %g\n",
		  init,ilay+1, glaya[ilay],N0_orig[ilay],
		  E->viscosity.N0[ilay],scale,E->viscosity.adjust_N0_alpha);
      }
    }
  }
  init++;
  return;
}


void read_visc_layer_file(struct All_variables *E)
{
    int i;
    FILE *fp;
    char junk[256];

    fp = fopen(E->viscosity.layer_file, "r");
    if (fp == NULL) {
        fprintf(E->fp, "(Viscosity_structures #1) Cannot open %s\n", E->viscosity.layer_file);
        exit(8);
    }


    /* default value */
    for(i=0; i<E->mesh.elz; i++) {
        E->viscosity.N0[i] =
            E->viscosity.E[i] =
            E->viscosity.T[i] =
  	    E->viscosity.Te[i] =
            E->viscosity.Z[i] =
            E->viscosity.pdepv_a[i] =
            E->viscosity.pdepv_b[i] =
            E->viscosity.pdepv_y[i] =
            E->viscosity.sdepv_expt[i] = 0;
    }

    read_visc_param_from_file(E, "visc0", E->viscosity.N0, fp);
    if(E->viscosity.TDEPV) {
        read_visc_param_from_file(E, "viscE", E->viscosity.E, fp);
        read_visc_param_from_file(E, "viscT", E->viscosity.T, fp);
	read_visc_param_from_file(E, "viscTe", E->viscosity.Te, fp);
        read_visc_param_from_file(E, "viscZ", E->viscosity.Z, fp);
    }

    if(E->viscosity.SDEPV) {
        read_visc_param_from_file(E, "sdepv_expt", E->viscosity.sdepv_expt, fp);
    }

    if(E->viscosity.PDEPV) {
        read_visc_param_from_file(E, "pdepv_a", E->viscosity.pdepv_a, fp);
        read_visc_param_from_file(E, "pdepv_b", E->viscosity.pdepv_b, fp);
        read_visc_param_from_file(E, "pdepv_y", E->viscosity.pdepv_y, fp);
    }

    return;
}


void visc_from_T(E,EEta,propogate)
     struct All_variables *E;
     float **EEta;
     int propogate;
{
    int m,i,k,l,z,jj,kk;
    float zero,one,eta0,temp,tempa,TT[9];
    float zzz,zz[9],dr;
    float visc1, visc2;
    const int vpts = vpoints[E->mesh.nsd];
    const int ends = enodes[E->mesh.nsd];
    const int nel = E->lmesh.nel;

    one = 1.0;
    zero = 0.0;

    /* consistent handling : l is (material number - 1) to allow
       addressing viscosity arrays, which are all 0...n-1  */
    switch (E->viscosity.RHEOL)   {
    case 1:
        /* eta = N_0 exp( E * (T_0 - T))  */
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;

                if(E->control.mat_control==0)
                    tempa = E->viscosity.N0[l];
                else 
                    tempa = E->viscosity.N0[l]*E->VIP[m][i];

                for(kk=1;kk<=ends;kk++) {
                    TT[kk] = E->T[m][E->ien[m][i].node[kk]];
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=0.0;
                    for(kk=1;kk<=ends;kk++)   {
                        temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)];
                    }

                    EEta[m][ (i-1)*vpts + jj ] = tempa*
                        exp( E->viscosity.E[l] * (E->viscosity.T[l] - temp));

                }
            }
        break;

    case 2:
        /* eta = N_0 exp(-T/T_0) */
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;

                if(E->control.mat_control==0)
                    tempa = E->viscosity.N0[l];
                else 
                    tempa = E->viscosity.N0[l]*E->VIP[m][i];

                for(kk=1;kk<=ends;kk++) {
                    TT[kk] = E->T[m][E->ien[m][i].node[kk]];
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=0.0;
                    for(kk=1;kk<=ends;kk++)   {
                        temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)];
                    }

                    EEta[m][ (i-1)*vpts + jj ] = tempa*
                        exp( -temp / E->viscosity.T[l]);

                }
            }
        break;

    case 3:
        /* eta = N_0 exp(E/(T+T_0) - E/(1+T_0)) 

	   where T is normalized to be within 0...1

	 */
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;
		if(E->control.mat_control) /* switch moved up here TWB */
		  tempa = E->viscosity.N0[l] * E->VIP[m][i];
		else
		  tempa = E->viscosity.N0[l];

                for(kk=1;kk<=ends;kk++) {
		  TT[kk] = E->T[m][E->ien[m][i].node[kk]];
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=0.0;
                    for(kk=1;kk<=ends;kk++)   {	/* took out
						   computation of
						   depth, not needed
						   TWB */
		      TT[kk]=max(TT[kk],zero);
		      temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
                    }
		    EEta[m][ (i-1)*vpts + jj ] = tempa*
		      exp( E->viscosity.E[l]/(temp+E->viscosity.T[l])
			   - E->viscosity.E[l]/(one +E->viscosity.T[l]) );
                }
            }
        break;

    case 4:
        /* eta = N_0 exp( (E + (1-z)Z_0) / (T+T_0) ) */
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;
		if(E->control.mat_control) /* moved this up here TWB */
		  tempa = E->viscosity.N0[l] * E->VIP[m][i];
		else
		  tempa = E->viscosity.N0[l];

                for(kk=1;kk<=ends;kk++) {
                    TT[kk] = E->T[m][E->ien[m][i].node[kk]];
                    zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=0.0;
                    zzz=0.0;
                    for(kk=1;kk<=ends;kk++)   {
                        TT[kk]=max(TT[kk],zero);
                        temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
                        zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
                    }


		    EEta[m][ (i-1)*vpts + jj ] = tempa*
		      exp( (E->viscosity.E[l] +  E->viscosity.Z[l]*zzz )
			   / (E->viscosity.T[l]+temp) );

                }
            }
        break;


    case 5:

        /* when mat_control=0, same as rheol 3,
           when mat_control=1, applying viscosity cut-off before mat_control */
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;
                tempa = E->viscosity.N0[l];
                /* fprintf(stderr,"\nINSIDE visc_from_T, l=%d, tempa=%g",l+1,tempa);*/
                for(kk=1;kk<=ends;kk++) {
                    TT[kk] = E->T[m][E->ien[m][i].node[kk]];
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=0.0;
                    for(kk=1;kk<=ends;kk++)   {
                        TT[kk]=max(TT[kk],zero);
                        temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
                    }

                    if(E->control.mat_control==0){
                        EEta[m][ (i-1)*vpts + jj ] = tempa*
                            exp( E->viscosity.E[l]/(temp+E->viscosity.T[l])
                                 - E->viscosity.E[l]/(one +E->viscosity.T[l]) );
		    }else{
                       visc2 = tempa*
	               exp( E->viscosity.E[l]/(temp+E->viscosity.T[l])
		          - E->viscosity.E[l]/(one +E->viscosity.T[l]) );
                       if(E->viscosity.MAX) {
                           if(visc2 > E->viscosity.max_value)
                               visc2 = E->viscosity.max_value;
                         }
                       if(E->viscosity.MIN) {
                           if(visc2 < E->viscosity.min_value)
                               visc2 = E->viscosity.min_value;
                         }
                       EEta[m][ (i-1)*vpts + jj ] = E->VIP[m][i]*visc2;
                      }

                }
            }
        break;


    case 6:
        /* like case 1, but allowing for depth-dependence if Z_0 != 0
           eta = N_0 exp(E(T_0-T) + (1-z) Z_0 )
        */

        for(m=1;m <= E->sphere.caps_per_proc;m++)
	  for(i=1;i <= nel;i++)   {

	    l = E->mat[m][i] - 1;

	    if(E->control.mat_control)
	      tempa = E->viscosity.N0[l] * E->VIP[m][i];
	    else
	      tempa = E->viscosity.N0[l];

	    for(kk=1;kk<=ends;kk++) {
	      TT[kk] = E->T[m][E->ien[m][i].node[kk]];
	      zz[kk] = (1.0 - E->sx[m][3][E->ien[m][i].node[kk]]);
	    }

	    for(jj=1;jj <= vpts;jj++) {
	      temp=0.0;zzz=0.0;
	      for(kk=1;kk <= ends;kk++)   {
		TT[kk]=max(TT[kk],zero);
		temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
		zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
	      }
	      EEta[m][ (i-1)*vpts + jj ] = tempa*
		exp( E->viscosity.E[l]*(E->viscosity.T[l] - temp) +
		     zzz *  E->viscosity.Z[l]);
	      
	      /* if(E->parallel.me == 0) */
	      /* 	fprintf(stderr,"z %11g km mat %i N0 %11g T %11g T0 %11g E %11g Z %11g eta: %12g\n", */
	      /* 		zzz *E->data.radius_km ,l+1, */
	      /* 		tempa,temp,E->viscosity.T[l],E->viscosity.E[l], E->viscosity.Z[l],EEta[m][ (i-1)*vpts + jj ]); */
	    }
	  }
        break;


    case 7:
        /* The viscosity formulation (dimensional) is:
           visc=visc0*exp[(Ea+p*Va)/(R*T)]

           Typical values for dry upper mantle are:
           Ea = 300 KJ/mol ; Va = 1.e-5 m^3/mol

           T=DT*(T0+T');
           where DT - temperature contrast (from Rayleigh number)
           T' - nondimensional temperature;
           T0 - nondimensional surface tempereture;

           =>
           visc = visc0 * exp{(Ea+p*Va) / [R*DT*(T0 + T')]}
                = visc0 * exp{[Ea/(R*DT) + p*Va/(R*DT)] / (T0 + T')}

           so:
           E->viscosity.E = Ea/(R*DT);
           (1-r) = p/(rho*g);
           E->viscosity.Z = Va*rho*g/(R*DT);
           E->viscosity.T = T0;

           after normalizing visc=1 at T'=1 and r=r_CMB:
           visc = visc0*exp{ [viscE + (1-r)*viscZ] / (viscT+T')
                - [viscE + (1-r_CMB)*viscZ] / (viscT+1) }
        */

        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
	      l = E->mat[m][i] - 1;

		if(E->control.mat_control)
		  tempa = E->viscosity.N0[l] * E->VIP[m][i];
		else
		  tempa = E->viscosity.N0[l];

                for(kk=1;kk<=ends;kk++) {
                    TT[kk] = E->T[m][E->ien[m][i].node[kk]];
                    zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=0.0;
                    zzz=0.0;
                    for(kk=1;kk<=ends;kk++)   {
                        temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)];
                        zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
                    }


                    EEta[m][ (i-1)*vpts + jj ] = tempa*
                        exp( (E->viscosity.E[l] +  E->viscosity.Z[l]*zzz )
                             / (E->viscosity.T[l] + temp)
                             - (E->viscosity.E[l] +
                                E->viscosity.Z[l]*(E->sphere.ro-E->sphere.ri) )
                             / (E->viscosity.T[l] + one) );
                }
            }
        break;

    case 8:
        /*
          eta0 = N_0 exp(E/(T+T_0) - E/(1+T_0))

          eta =        eta0 if T  < T_sol0 + 2(1-z)
          eta = ET_red*eta0 if T >= T_sol0 + 2(1-z)

	  T is limited to lie between 0 and 1

          where z is normalized by layer
          thickness, and T_sol0 is something
          like 0.6, and ET_red = 0.1

          (same as case 3, but for viscosity reduction)
        */

        dr = E->sphere.ro - E->sphere.ri;
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;
		if(E->control.mat_control) 
		  tempa = E->viscosity.N0[l] * E->VIP[m][i];
		else
		  tempa = E->viscosity.N0[l];

                for(kk=1;kk<=ends;kk++) {
		  TT[kk] = E->T[m][E->ien[m][i].node[kk]];
		  zz[kk] = E->sx[m][3][E->ien[m][i].node[kk]]; /* radius */
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=zzz=0.0;
                    for(kk=1;kk<=ends;kk++)   {	
		      TT[kk]=max(TT[kk],zero);
		      temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)]; /* mean temp */
		      zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];/* mean r */
                    }
		    /* convert to z, as defined to be unity at surface
		       and zero at CMB */
		    zzz = (zzz - E->sphere.ri)/dr;
		    visc1 = tempa* exp( E->viscosity.E[l]/(temp+E->viscosity.T[l]) 
				  - E->viscosity.E[l]/(one +E->viscosity.T[l]) );
		    if(temp < E->viscosity.T_sol0 + 2.*(1.-zzz))
		      EEta[m][ (i-1)*vpts + jj ] = visc1;
		    else
		      EEta[m][ (i-1)*vpts + jj ] = visc1 * E->viscosity.ET_red;
                }
            }
        break;
    case 9:
        /* eta = N_0 exp(E/(T+T_0) - E/(1+T_0)) 

	   like option 3, but T is allow to vary beyond 1 

	 */
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;
		if(E->control.mat_control) /* switch moved up here TWB */
		  tempa = E->viscosity.N0[l] * E->VIP[m][i];
		else
		  tempa = E->viscosity.N0[l];
                for(kk=1;kk<=ends;kk++) 
		  TT[kk] = E->T[m][E->ien[m][i].node[kk]];
		
                for(jj=1;jj<=vpts;jj++) {
                    temp=0.0;
                    for(kk=1;kk<=ends;kk++)
		      temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)];
		    EEta[m][ (i-1)*vpts + jj ] = tempa*
		      exp( E->viscosity.E[l]/(temp+E->viscosity.T[l])
			   - E->viscosity.E[l]/(one +E->viscosity.T[l]) );
                }
            }
        break;
    case 10:
        /*
          eta0 = N_0 exp(E/(T+T_0) - E/(1+T_0))

          eta =        eta0 if T  < T_sol0 + 2(1-z)
          eta = ET_red*eta0 if T >= T_sol0 + 2(1-z)

	  like rheol == 8, but T is not limited to lie between 0 and 1

        */

        dr = E->sphere.ro - E->sphere.ri;
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;
		if(E->control.mat_control) 
		  tempa = E->viscosity.N0[l] * E->VIP[m][i];
		else
		  tempa = E->viscosity.N0[l];

                for(kk=1;kk<=ends;kk++) {
		  TT[kk] = E->T[m][E->ien[m][i].node[kk]];
		  zz[kk] = E->sx[m][3][E->ien[m][i].node[kk]]; /* radius */
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=zzz=0.0;
                    for(kk=1;kk<=ends;kk++)   {	
		      temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)]; /* mean temp */
		      zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];/* mean r */
                    }
		    zzz = (zzz - E->sphere.ri)/dr;
		    visc1 = tempa* exp( E->viscosity.E[l]/(temp+E->viscosity.T[l]) 
				  - E->viscosity.E[l]/(one +E->viscosity.T[l]) );
		    if(temp < E->viscosity.T_sol0 + 2.*(1.-zzz))
		      EEta[m][ (i-1)*vpts + jj ] = visc1;
		    else
		      EEta[m][ (i-1)*vpts + jj ] = visc1 * E->viscosity.ET_red;
                }
            }
        break;

    case 11:
        /* eta = N_0 exp(E/(T+T_0) - E/(T_eta+T_0)) 
 *		
 *		where T_0 is the offset temperature
 *		and T_eta the reference temperature for eta(T) = 1
 *		where T is normalized to be within 0...1
 *
 *                     */
        for(m=1;m<=E->sphere.caps_per_proc;m++)
            for(i=1;i<=nel;i++)   {
                l = E->mat[m][i] - 1;
                if(E->control.mat_control) /* switch moved up here TWB */
                  tempa = E->viscosity.N0[l] * E->VIP[m][i];
                else
                  tempa = E->viscosity.N0[l];

                for(kk=1;kk<=ends;kk++) {
                  TT[kk] = E->T[m][E->ien[m][i].node[kk]];
                }

                for(jj=1;jj<=vpts;jj++) {
                    temp=0.0;
                    for(kk=1;kk<=ends;kk++)   { /* took out
                                                   computation of
                                                   depth, not needed
                                                   TWB */
                      TT[kk]=max(TT[kk],zero);
                      temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
                    }
                    EEta[m][ (i-1)*vpts + jj ] = tempa*
                      exp( E->viscosity.E[l]/(temp+E->viscosity.T[l])
                           - E->viscosity.E[l]/(E->viscosity.Te[l]+E->viscosity.T[l]) );
                }
            }
        break;	

    case 100:
        /* user-defined viscosity law goes here */
        fprintf(stderr, "Need user definition for viscosity law: 'rheol=%d'\n",
                E->viscosity.RHEOL);
        parallel_process_termination();
        break;

    default:
        /* unknown option */
        fprintf(stderr, "Invalid value of 'rheol=%d'\n", E->viscosity.RHEOL);

        parallel_process_termination();
        break;
    }


    return;
}
/* compute the horizontally averaged strain-rates, edot[nelz+1] */
void calc_radial_strain_rate_averages(struct All_variables *E, float *edot)
{
   const int nel = E->lmesh.nel;
   int m;
   float *eedot;
   eedot = (float *) malloc((2+nel)*sizeof(float));
   for(m=1;m<=E->sphere.caps_per_proc;m++)  {
     strain_rate_2_inv(E,m,eedot,1);
     return_elementwise_horiz_ave_float(E,m,eedot,edot);
   }
   free(eedot);
}

void visc_from_S(E,EEta,propogate)
     struct All_variables *E;
     float **EEta;
     int propogate;
{
  float scale,exponent1;
  float *eedot;
  int m,e,l,z,jj,kk;

  const int vpts = vpoints[E->mesh.nsd];
  const int nel = E->lmesh.nel;
  
  eedot = (float *) malloc((2+nel)*sizeof(float));
  
  for(m=1;m<=E->sphere.caps_per_proc;m++)  {
    if(E->viscosity.sdepv_visited){
      
      /* get second invariant for all elements */
      strain_rate_2_inv(E,m,eedot,1);
    }else{
      for(e=1;e<=nel;e++)	/* initialize with unity if no velocities around */
	eedot[e] = 1.0;
      E->viscosity.sdepv_visited = 1;
      
    }
    /* eedot cannot be too small, or the viscosity will go to inf */
    for(e=1;e<=nel;e++){
      eedot[e] = max(eedot[e], 1.0e-16);
    }
    if(!E->viscosity.sdepv_start_from_newtonian){
      switch(E->viscosity.sdepv_rheol){
      case 1:		/* old default - not sure this is any good */ 
	for(e=1;e<=nel;e++)   {
	  exponent1= 1.0/E->viscosity.sdepv_expt[E->mat[m][e]-1];
	  scale=pow(eedot[e],exponent1-1.0);
	  for(jj=1;jj<=vpts;jj++)
	    EEta[m][(e-1)*vpts + jj] = scale*pow(EEta[m][(e-1)*vpts+jj],exponent1);
	}
	break;
      case 2:			/* new */
	for(e=1;e<=nel;e++)   {
	  exponent1 = 1.0 - 1.0 / E->viscosity.sdepv_expt[E->mat[m][e] - 1];
	  scale = pow(2.0 * eedot[e] / E->viscosity.sdepv_trns[E->mat[m][e] - 1], exponent1);
	  //if(E->parallel.me==0){      /* control output */
	  // fprintf(stderr,"Power-law, strain = %e, trans_strain = %e, exp1 = %e  ET = %e ETP = %e \n",
	  //	    eedot[e],E->viscosity.sdepv_trns[E->mat[m][e] - 1],exponent1,EEta[m][(e - 1) * vpts + 1],
	  //	    EEta[m][(e - 1) * vpts + 1] / (1.0 + scale * pow(EEta[m][(e - 1) * vpts + 1], exponent1)));
	  //}
	  for(jj = 1; jj <= vpts; jj++){
	    EEta[m][(e - 1) * vpts + jj] = 
	      EEta[m][(e - 1) * vpts + jj] / (1.0 + scale * pow(EEta[m][(e - 1) * vpts + jj], exponent1));
	  }
	}
	break;
      }
      /* get second invariant for all elements */
      strain_rate_2_inv(E,m,eedot,1);
    }else{
      for(e=1;e<=nel;e++)	/* initialize with unity if no velocities around */
	eedot[e] = 1.0;
      E->viscosity.sdepv_visited = 1;
      
    }
   
  }
  free ((void *)eedot);
  return;
}

void visc_from_P(E,EEta) /* "plasticity" implementation


			 psrw = FALSE

			 viscosity will be limited by a yield stress
			 
			 \sigma_y  = min(a + b * (1-r), y)

			 where a,b,y are parameters input via pdepv_a,b,y

			 and

			 \eta_y = \sigma_y / (2 \eps_II)

			 where \eps_II is the second invariant. Then

			 \eta_eff = (\eta_0 \eta_y)/(\eta_0 + \eta_y)

			 for pdepv_eff = 1

			 or

			 \eta_eff = min(\eta_0,\eta_y)

			 for pdepv_eff = 0

			 where \eta_0 is the regular viscosity


			 psrw = 1

			 a strain-rate weakening rheology applies
			 based on a steady state stress-strain
			 relationship following, e.g., Tackley 1998

			 \eta_ett = (\sigma_y^2 \eta_0)/(\sigma_y^2 + \eta_0^2 (2 \eps_II^2))

			 where sigma_y is defined as above

			 where the 2\eps_II arises because our \eps_II has the 1/2 factor in it

			 TWB

			 */
     struct All_variables *E;
     float **EEta;
{
  float *eedot,zz[9],zzz,tau,eta_p,eta_new,tau2,eta_old,eta_old2,ss[9],sss;
  float fac;
  double estrain,maxD_over_gamma_cr,peak_strain_rate,scaled_strain_rate;
  int m,e,l,z,jj,kk;
  float *eetty[NCS], *esigY[NCS]; 
  int check_stress_state;
  const int vpts = vpoints[E->mesh.nsd];
  const int nel = E->lmesh.nel;
  const int ends = enodes[E->mesh.nsd];
  
  void strain_rate_2_inv();

  maxD_over_gamma_cr = E->viscosity.maxD/E->viscosity.gamma_cr;
  
  if(E->viscosity.p_stress_state_dependent && E->computed_gstress){
    check_stress_state = E->viscosity.p_stress_state_dependent; /* 1: for weakening 2: for plasticity */
  }else{
    check_stress_state = 0;
  }
  
  eedot = (float *) malloc((2+nel)*sizeof(float));

  for(m=1;m<=E->sphere.caps_per_proc;m++) {
    eetty[m] = (float *) malloc((2+nel)*sizeof(float));
    esigY[m] = (float *) malloc((2+nel)*sizeof(float));
  }

  if(E->control.restart) {
     E->viscosity.pdepv_visited = 1;
  }

  for(m=1;m<=E->sphere.caps_per_proc;m++) {
    
    if(E->viscosity.pdepv_visited){
      if(E->viscosity.psrw)
	strain_rate_2_inv(E,m,eedot,0);	/* get second invariant for all elements, don't take sqrt */
      else
	strain_rate_2_inv(E,m,eedot,1);	/* get second invariant for all elements */
    }else{
      for(e=1;e<=nel;e++)	/* initialize with unity if no velocities around */
	eedot[e] = 1.0;
      if(m == E->sphere.caps_per_proc)
	E->viscosity.pdepv_visited = 1;
      if((E->parallel.me == 0)&&(E->control.verbose)){
	fprintf(stderr,"num mat: %i a: %g b: %g y: %g %s\n",
		e,E->viscosity.pdepv_a[e],E->viscosity.pdepv_b[e],E->viscosity.pdepv_y[e],
		(E->viscosity.psrw)?(" -- SRW"):(""));
      }
    }
    if(!E->viscosity.psrw){
      /* 
	 regular plasticity
      */
      for(e=1;e <= nel;e++)   {	/* loop through all elements */
	
	l = E->mat[m][e] -1 ;	/* material of this element */
	
	for(kk=1;kk <= ends;kk++) {
	  /* nodal depths */
	  zz[kk] = (1.0 - E->sx[m][3][E->ien[m][e].node[kk]]); /* for depth, zz = 1 - r */
	  /* how strike slip is this node? 1: pure 0: normal/extensional */ 
	  if(check_stress_state){
	    ss[kk] = compute_nodal_stress_state_factor(E,m,E->ien[m][e].node[kk]); 
	  }
	}
        eetty[m][e]=0.0;
        esigY[m][e]=0.0; 

	for(jj=1;jj <= vpts;jj++){ /* loop through integration points */
	  
	  zzz = 0.0;		/* get mean depth of integration point */
	  if(check_stress_state)
	    sss = 0.0;
	  else
	    sss = 1.0;		/* default, use full damage */
	  for(kk=1;kk<=ends;kk++){
	    zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
	    if(check_stress_state){
	      sss += ss[kk] * E->N.vpt[GNVINDEX(kk,jj)]; /* state, if > 0.5, will weaken */
	    }
	  }
	  
	  /* depth dependent yield stress */
	  tau = E->viscosity.pdepv_a[l] + zzz * E->viscosity.pdepv_b[l];
	  
	  /* min of depth dep. and constant yield stress */
	  tau = min(tau,  E->viscosity.pdepv_y[l]);

	  if(check_stress_state == 2){
	    /* 
	       apply stress state to regular plasticity . 2....0.5

	     */
	    tau *= pow(10,tanh((0.5-sss)*7)/5);
	    sss = 1;		/* make sure all the other weakening
				   applies */
	  }


          /* Add compositional contribution? */
          if(E->control.stiff_cont){
             /* Tau = tau for C = 0, tau = tau_c for C = 1, and linearly in
                between  */
             tau *= (1.0 - E->composition.comp_el[m][0][e]);
             tau += E->composition.comp_el[m][0][e] * E->viscosity.sig_c;
	  }

	  /* 
	     Add strain dependent plasticity function 
	  */
          if(E->viscosity.strain_dep_plasticity == 1){
	    // "Get" strain from element
	    estrain = E->trace.strain_el[m][e] * sss;
	    // Weakening term 
	    fac = 1.0 - estrain * maxD_over_gamma_cr;
	    if(fac < (1.0-E->viscosity.maxD))
	      fac = (1.0-E->viscosity.maxD);
	    // Reduced yield stress
	    tau *= fac;
	  }

	  esigY[m][e]=tau;

	  /* yield viscosity */
	  eta_p = tau/(2.0 * eedot[e] + 1e-7) + E->viscosity.pdepv_offset;

	  /* calculate and store yield viscosity ratio 
		y_fac = eta(T)/eta_p
          */

          if(EEta[m][ (e-1)*vpts + jj ] > E->viscosity.max_value){
            eetty[m][e] += E->viscosity.max_value/eta_p;
          }else if(EEta[m][ (e-1)*vpts + jj ] < E->viscosity.min_value){
            eetty[m][e] += E->viscosity.min_value/eta_p;
          }else{
            eetty[m][e] += EEta[m][ (e-1)*vpts + jj ]/eta_p;
          }
          
          if(E->viscosity.strain_dep_plasticity == 2){
	    estrain = E->trace.strain_el[m][e] * sss;
	    // Weakening term
	    fac = 1.0 - estrain * maxD_over_gamma_cr;
	    if(fac < (1.0-E->viscosity.maxD))
	      fac = (1.0-E->viscosity.maxD);
	    EEta[m][ (e-1)*vpts + jj ] *= fac;
	  }

	  if(E->viscosity.pdepv_eff){
	    /* two dashpots in series */
	    eta_new  = 1.0/(1.0/EEta[m][ (e-1)*vpts + jj ] + 1.0/eta_p);
	  }else{
	    /* min viscosities*/
	    eta_new  = min(EEta[m][ (e-1)*vpts + jj ], eta_p);
	  }

          if(E->viscosity.strain_dep_plasticity == 3){
	    estrain = E->trace.strain_el[m][e] * sss;
	    // Weakening term
	    fac = 1.0 - estrain * maxD_over_gamma_cr;
	    if(fac < (1.0-E->viscosity.maxD))
	      fac = (1.0-E->viscosity.maxD);
	    // eta_eff = eta_eff * (1-D)^Delta_etaT; 
	    eta_new *= pow(fac,exp( E->viscosity.E[l]/2.0));
	  }
          
	  //fprintf(stderr,"z: %11g mat: %i a: %11g b: %11g y: %11g ee: %11g tau: %11g eta_p: %11g eta_new: %11g eta_old: %11g\n",
	  //	  zzz,l,E->viscosity.pdepv_a[l], E->viscosity.pdepv_b[l],E->viscosity.pdepv_y[l],
	  //	  eedot[e],tau,eta_p,eta_new,EEta[m][(e-1)*vpts + jj]);
	  EEta[m][(e-1)*vpts + jj] = eta_new;
	} /* end integration point loop */
        eetty[m][e]/=vpts; 
      }	/* end element loop */
    }else{
      /* strain-rate weakening, steady state solution 

	 \tau = 2/e_c (e/(1+(e/e_c)^2))
	 
	 where e_c is the strain-rate at maximum stress, here \sigma_y, hence e_c = \sigma_y/(2\eta_T)
	 
       */
      for(e=1;e <= nel;e++)   {	/* loop through all elements */
	
	l = E->mat[m][e] -1 ;	
	for(kk=1;kk <= ends;kk++)
	  zz[kk] = (1.0 - E->sx[m][3][E->ien[m][e].node[kk]]); 
	for(jj=1;jj <= vpts;jj++){ 
	  zzz = 0.0;
	  for(kk=1;kk<=ends;kk++)
	    zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
	  /* compute sigma_y as above */
	  tau = E->viscosity.pdepv_a[l] + zzz * E->viscosity.pdepv_b[l];
	  tau = min(tau,  E->viscosity.pdepv_y[l]);
	  /*  */
	  eta_old = EEta[m][ (e-1)*vpts + jj ];
	  /*  */
	  peak_strain_rate = tau / (2.0*eta_old+1e-15);
	  /* effective viscosity */
	  scaled_strain_rate = eedot[e]/peak_strain_rate;
	  
	  eta_new = eta_old / (peak_strain_rate * (1.0 +  scaled_strain_rate *  scaled_strain_rate));

	  //fprintf(stderr,"SRW: a %11g b %11g y %11g z %11g sy: %11g e2: %11g eold: %11g enew: %11g logr: %.3f\n",
	  //	    E->viscosity.pdepv_a[l],E->viscosity.pdepv_b[l],E->viscosity.pdepv_y[l],zzz,tau,eedot[e],eta_old,eta_new,
	  //	    log10(eta_new/eta_old));
	  EEta[m][(e-1)*vpts + jj] = min(eta_new, eta_old);
	}
      }
    }
  } /* end caps loop */
  /* project yield viscosity from element to nodes */
  viscy_from_element_to_nodes(E,eetty,E->VIY[E->mesh.levmax],E->mesh.levmax);
  viscy_from_element_to_nodes(E,esigY,E->SIGY[E->mesh.levmax],E->mesh.levmax);
  
  free ((void *)eedot);

  for(m=1;m<=E->sphere.caps_per_proc;m++){
    free ((void *)eetty[m]); 
    free ((void *)esigY[m]); 
  }
  return;
}

/*

multiply with compositional factor which is determined by a geometric
mean average from the tracer composition, assuming two flavors and
compositions between zero and unity

*/
void visc_from_C( E, EEta)
     struct All_variables *E;
     float **EEta;
{
  double vmean,cc_loc[10],CC[10][9],cbackground;
  int m,l,z,jj,kk,i,p,q;


  const int vpts = vpoints[E->mesh.nsd];
  const int nel = E->lmesh.nel;
  const int ends = enodes[E->mesh.nsd];

  for(m=1;m <= E->sphere.caps_per_proc;m++)  {
    for(i = 1; i <= nel; i++){
      /* determine composition of each of the nodes of the
	 element */
        for(p=0; p<E->composition.ncomp; p++) {
            for(kk = 1; kk <= ends; kk++){
                CC[p][kk] = E->composition.comp_node[m][p][E->ien[m][i].node[kk]];
                if(CC[p][kk] < 0)CC[p][kk]=0.0;
                if(CC[p][kk] > 1)CC[p][kk]=1.0;
            }
        }
        for(jj = 1; jj <= vpts; jj++) {
            /* concentration of background material */
            cbackground = 1;
            for(p=0; p<E->composition.ncomp; p++) {
                /* compute mean composition  */
                cc_loc[p] = 0.0;
                for(kk = 1; kk <= ends; kk++) {
                    cc_loc[p] += CC[p][kk] * E->N.vpt[GNVINDEX(kk, jj)];
                }
                cbackground -= cc_loc[p];
            }

            /* geometric mean of viscosity */
            vmean = cbackground * E->viscosity.cdepv_ff[0];
            for(p=0; p<E->composition.ncomp; p++) {
                vmean += cc_loc[p] * E->viscosity.cdepv_ff[p+1];
            }
            vmean = exp(vmean);

            /* multiply the viscosity with this prefactor */
            EEta[m][ (i-1)*vpts + jj ] *= vmean;

        } /* end jj loop */
    } /* end el loop */
  } /* end cap */
}

/* 

   compute the stress state that determines how much of the damage applies

*/
float compute_nodal_stress_state_factor(struct All_variables *E, int m, int node)
{
  float a[3][3],tnorm,state,meanh;
  int i,j;
  a[0][0]=E->gstress[m][(node-1)*6+3]; /*  srr */
  a[1][1]=E->gstress[m][(node-1)*6+1]; /*  stt */
  a[2][2]=E->gstress[m][(node-1)*6+2]; /*  spp */
  a[0][1]=a[1][0]=E->gstress[m][(node-1)*6+5]; /*  str */
  a[0][2]=a[2][0]=E->gstress[m][(node-1)*6+6]; /* srp */
  a[1][2]=a[2][1]=E->gstress[m][(node-1)*6+4]; /*  stp */


  tnorm = 0;				       /* tensor norm */
  for(i=0;i<3;i++)
    for(j=0;j<3;j++)
      tnorm += a[i][j] * a[i][j];
  tnorm = sqrt(tnorm);
  meanh = (a[1][1]+a[2][2])/2/tnorm; /* this goes from ~-0.71 = sqrt(2)/2
					(compressional) to 0.71
					(extensional) in the
					horizontal */
  state = 1-fabs(meanh)/0.35;
  if(state<0)
    state = 0;
  return state;			/* 1=strike slip 0: compressional or extensional */
  
}

void strain_rate_2_inv(E,m,EEDOT,SQRT)
     struct All_variables *E;
     float *EEDOT;
     int m,SQRT;
{
    void get_rtf_at_ppts();
    void velo_from_element();
    void construct_c3x3matrix_el();
    void get_ba_p();

    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;

    const int nel = E->lmesh.nel;
    const int dims = E->mesh.nsd;
    const int ends = enodes[dims];
    const int lev = E->mesh.levmax;
    const int ppts = ppoints[dims];
    const int sphere_key = 1;

    for(e=1; e<=nel; e++) {

        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]);

        theta = rtf[1][1];


        /* Vxyz is the strain rate vector, whose relationship with
         * the strain rate tensor (e) is that:
         *    Vxyz[1] = e11
         *    Vxyz[2] = e22
         *    Vxyz[3] = e33
         *    Vxyz[4] = 2*e12
         *    Vxyz[5] = 2*e13
         *    Vxyz[6] = 2*e23
         * where 1 is theta, 2 is phi, and 3 is r
         */
        for(j=1; j<=ppts; j++) {
            Vxyz[1][j] = 0.0;
            Vxyz[2][j] = 0.0;
            Vxyz[3][j] = 0.0;
            Vxyz[4][j] = 0.0;
            Vxyz[5][j] = 0.0;
            Vxyz[6][j] = 0.0;
            dilation[j] = 0.0;
        }

        if ((E->control.precise_strain_rate) || (theta < 0.09) || (theta > 3.05)) {
            /* When the element is close to the poles, use a more
             * precise method to compute the strain rate. 
	     
	     if precise_strain_rate=on, will always choose this option

	    */

            if ((e-1)%E->lmesh.elz==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][q][p] * 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 fast, imprecise strain-rate computation */

        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;
        }

        edot[1][1] = edot[2][2] = edot[3][3] = 0;
        edot[1][2] = edot[1][3] = edot[2][3] = 0;

        /* edot is 2 * (the deviatoric strain rate tensor) */
        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]);
            edot[1][2] += Vxyz[4][j];
            edot[1][3] += Vxyz[5][j];
            edot[2][3] += Vxyz[6][j];
        }

        EEDOT[e] = edot[1][1] * edot[1][1]
            + edot[1][2] * edot[1][2] * 2.0
            + edot[2][2] * edot[2][2]
            + edot[2][3] * edot[2][3] * 2.0
            + edot[3][3] * edot[3][3]
            + edot[1][3] * edot[1][3] * 2.0;
    }

    if(SQRT)
	for(e=1;e<=nel;e++)
	    EEDOT[e] =  sqrt(0.5 *EEDOT[e]);
    else
	for(e=1;e<=nel;e++)
	    EEDOT[e] *=  0.5;

    return;
}


static void apply_low_visc_wedge_channel(struct All_variables *E, float **evisc)
{
    void parallel_process_termination();

    int i,j,m;
    const int vpts = vpoints[E->mesh.nsd];
    float *F;

    /* low viscosity channel/wedge require tracers to work */
    if(E->control.tracer == 0) {
        if(E->parallel.me == 0) {
            fprintf(stderr, "Error: low viscosity channel/wedge is turned on, "
                   "but tracer is off!\n");
            fprintf(E->fp, "Error: low viscosity channel/wedge is turned on, "
                   "but tracer is off!\n");
            fflush(E->fp);
        }
        parallel_process_termination();
    }


    F = (float *)malloc((E->lmesh.nel+1)*sizeof(float));
    for(i=1 ; i<=E->lmesh.nel ; i++)
        F[i] = 0.0;

    /* if low viscosity channel ... */
    if(E->viscosity.channel)
        low_viscosity_channel_factor(E, F);


    /* if low viscosity wedge ... */
    if(E->viscosity.wedge)
        low_viscosity_wedge_factor(E, F);


    for(i=1 ; i<=E->lmesh.nel ; i++) {
        if (F[i] != 0.0)
            for(m = 1 ; m <= E->sphere.caps_per_proc ; m++) {
                for(j=1;j<=vpts;j++) {
                    evisc[m][(i-1)*vpts + j] = F[i];
            }
        }
    }


    free(F);

    return;
}




static void low_viscosity_channel_factor(struct All_variables *E, float *F)
{
    int i, ii, k, m, e, ee;
    int nz_min[NCS], nz_max[NCS];
    const int flavor = 0;
    double rad_mean, rr;

    for(m=1; m<=E->sphere.caps_per_proc; m++) {
        /* find index of radius corresponding to lv_min_radius */
        for(e=1; e<=E->lmesh.elz; e++) {
            rad_mean = 0.5 * (E->sx[m][3][E->ien[m][e].node[1]] +
                              E->sx[m][3][E->ien[m][e].node[8]]);
            if(rad_mean >= E->viscosity.lv_min_radius) break;
        }
        nz_min[m] = e;

        /* find index of radius corresponding to lv_max_radius */
        for(e=E->lmesh.elz; e>=1; e--) {
            rad_mean = 0.5 * (E->sx[m][3][E->ien[m][e].node[1]] +
                              E->sx[m][3][E->ien[m][e].node[8]]);
            if(rad_mean <= E->viscosity.lv_max_radius) break;
        }
        nz_max[m] = e;
    }



    for(m=1; m<=E->sphere.caps_per_proc; m++) {
        for(k=1; k<=E->lmesh.elx*E->lmesh.ely; k++) {
            for(i=nz_min[m]; i<=nz_max[m]; i++) {
                e = (k-1)*E->lmesh.elz + i;

                rad_mean = 0.5 * (E->sx[m][3][E->ien[m][e].node[1]] +
                                  E->sx[m][3][E->ien[m][e].node[8]]);

                /* loop over elements below e */
                for(ii=i; ii>=nz_min[m]; ii--) {
                    ee = (k-1)*E->lmesh.elz + ii;

                    rr = 0.5 * (E->sx[m][3][E->ien[m][ee].node[1]] +
                                E->sx[m][3][E->ien[m][ee].node[8]]);

                    /* if ee has tracers in it and is within the channel */
                    if((E->trace.ntracer_flavor[m][flavor][ee] > 0) &&
                       (rad_mean <= rr + E->viscosity.lv_channel_thickness)) {
                           F[e] = E->viscosity.lv_reduction;
                           break;
                       }
                }
            }
        }
    }


    return;
}


static void low_viscosity_wedge_factor(struct All_variables *E, float *F)
{
    int i, ii, k, m, e, ee;
    int nz_min[NCS], nz_max[NCS];
    const int flavor = 0;
    double rad_mean, rr;

    for(m=1; m<=E->sphere.caps_per_proc; m++) {
        /* find index of radius corresponding to lv_min_radius */
        for(e=1; e<=E->lmesh.elz; e++) {
            rad_mean = 0.5 * (E->sx[m][3][E->ien[m][e].node[1]] +
                              E->sx[m][3][E->ien[m][e].node[8]]);
            if(rad_mean >= E->viscosity.lv_min_radius) break;
        }
        nz_min[m] = e;

        /* find index of radius corresponding to lv_max_radius */
        for(e=E->lmesh.elz; e>=1; e--) {
            rad_mean = 0.5 * (E->sx[m][3][E->ien[m][e].node[1]] +
                              E->sx[m][3][E->ien[m][e].node[8]]);
            if(rad_mean <= E->viscosity.lv_max_radius) break;
        }
        nz_max[m] = e;
    }



    for(m=1; m<=E->sphere.caps_per_proc; m++) {
        for(k=1; k<=E->lmesh.elx*E->lmesh.ely; k++) {
            for(i=nz_min[m]; i<=nz_max[m]; i++) {
                e = (k-1)*E->lmesh.elz + i;

                rad_mean = 0.5 * (E->sx[m][3][E->ien[m][e].node[1]] +
                                  E->sx[m][3][E->ien[m][e].node[8]]);

                /* loop over elements below e */
                for(ii=i; ii>=nz_min[m]; ii--) {
                    ee = (k-1)*E->lmesh.elz + ii;

                    /* if ee has tracers in it */
                    if(E->trace.ntracer_flavor[m][flavor][ee] > 0) {
                        F[e] = E->viscosity.lv_reduction;
                        break;
                    }
                }
            }
        }
    }


    /* * debug **

    for(m=1; m<=E->sphere.caps_per_proc; m++)
        for(e=1; e<=E->lmesh.nel; e++)
            fprintf(stderr, "lv_reduction: %d %e\n", e, F[e]);
    */

    return;
}
/* compute second invariant from a strain-rate tensor in 0,...2 format

 */
double second_invariant_from_3x3(double e[3][3])
{
  return(sqrt(0.5*
	      (e[0][0] * e[0][0] + 
	       e[0][1] * e[0][1] * 2.0 + 
	       e[1][1] * e[1][1] + 
	       e[1][2] * e[1][2] * 2.0 + 
	       e[2][2] * e[2][2] + 
	       e[0][2] * e[0][2] * 2.0)));
}
void calc_strain_from_vgm(double l[3][3], double d[3][3])
{
  int i,j;
  for(i=0;i < 3;i++)
    for(j=0;j < 3;j++)
      d[i][j] = 0.5 * (l[i][j] + l[j][i]);
}
void calc_strain_from_vgm9(double *l9, double d[3][3])
{
  double l[3][3];
  get_3x3_from_9vec(l, l9);
  calc_strain_from_vgm(l, d);
}

/* 

   given a 3x3 velocity gradient matrix l, compute a rotation matrix

*/

void calc_rot_from_vgm(double l[3][3], double r[3][3])
{
  int i,j;
  for(i=0;i < 3;i++)
    for(j=0;j < 3;j++)
      r[i][j] = 0.5 * (l[i][j] - l[j][i]);
}

/* 

   get velocity gradient matrix at element, and also compute the
   average velocity in this element
   

*/

void get_vgm_p(double VV[4][9],struct Shape_function *N,
	       struct Shape_function_dx *GNx,
	       struct CC *cc, struct CCX *ccx, double rtf[4][9],
	       int dims,int ppts, int ends, int spherical,
	       double l[3][3], double v[3])
{

  int i,k,j,a;
  double ra[9], si[9], ct[9];
  const double one = 1.0;
  const double two = 2.0;
  double vgm[3][3];
  double shp, cc1, cc2, cc3,d_t,d_r,d_p,up,ur,ut;
  /* init L matrix */
  for(i=0;i < 3;i++){
    v[i] = 0.0;
    for(j=0;j < 3; j++)
      l[i][j] = 0.0;
  }
  /* mean velocity at pressure integration point */
  for(a=1;a <= ends;a++){
    v[0] += N->ppt[GNPINDEX(a, 1)] * VV[1][a];
    v[1] += N->ppt[GNPINDEX(a, 1)] * VV[2][a];
    v[2] += N->ppt[GNPINDEX(a, 1)] * VV[3][a];
  }
  if(spherical){
    for(k = 1; k <= ppts; k++){
      ra[k] = rtf[3][k];	      /* 1/r */
      si[k] = one / sin(rtf[1][k]); /* 1/sin(t) */
      ct[k] = cos(rtf[1][k]) * si[k]; /* cot(t) */
    }
    for(a = 1; a <= ends; a++){
      for(k = 1; k <= ppts; k++){
	d_t = GNx->ppt[GNPXINDEX(0, a, k)]; /* d_t */
	d_p = GNx->ppt[GNPXINDEX(1, a, k)]; /* d_p */
	d_r = GNx->ppt[GNPXINDEX(2, a, k)]; /* d_r */
	shp = N->ppt[GNPINDEX(a, k)];
	for(i = 1; i <= dims; i++){
	  ut = cc->ppt[BPINDEX(1, i, a, k)]; /* ut */
	  up = cc->ppt[BPINDEX(2, i, a, k)]; /* up */
	  ur = cc->ppt[BPINDEX(3, i, a, k)]; /* ur */
	  
	  /* velocity gradient matrix is transpose of grad v, using Citcom sort t, p, r
	
	     | d_t(vt) d_p(vt) d_r(vt) |
	     | d_t(vp) d_p(vp) d_r(vp) |
	     | d_t(vr) d_p(vr) d_r(vr) |

	  */

	  /* d_t vt = 1/r (d_t vt + vr) */
	  vgm[0][0] =  ((d_t * ut + shp * ccx->ppt[BPXINDEX(1, i, 1, a, k)]) + 
			shp * ur) * ra[k];
	  /* d_p vt = 1/r (1/sin(t) d_p vt -vp/tan(t)) */
	  vgm[0][1] =  ((d_p * ut + shp * ccx->ppt[BPXINDEX(1, i, 2, a, k)]) * si[k] - 
			shp * up * ct[k]) * ra[k];
	  /* d_r vt = d_r v_t */
	  vgm[0][2] = d_r * ut;
	  /* d_t vp = 1/r d_t v_p*/
	  vgm[1][0] = (d_t * up + shp * ccx->ppt[BPXINDEX(2, i, 1, a, k)]) * ra[k];
	  /* d_p vp = 1/r((d_p vp)/sin(t) + vt/tan(t) + vr) */
	  vgm[1][1] = ((d_p * up + shp * ccx->ppt[BPXINDEX(2, i, 2, a, k)]) * si[k] + 
		       shp * ut * ct[k] + shp * ur) * ra[k];
	  /* d_r vp = d_r v_p */
	  vgm[1][2] =  d_r * up;
	  /* d_t vr = 1/r(d_t vr - vt) */
	  vgm[2][0] = ((d_t * ur + shp * ccx->ppt[BPXINDEX(3, i, 1, a, k)]) -
		       shp * ut) * ra[k];
	  /* d_p vr =  1/r(1/sin(t) d_p vr - vp) */
	  vgm[2][1] = (( d_p * ur + shp * ccx->ppt[BPXINDEX(3,i, 2,a,k)] ) * si[k] -
		       shp * up ) * ra[k];
	  /* d_r vr = d_r vr */
	  vgm[2][2] = d_r * ur;


	  l[0][0] += vgm[0][0] * VV[i][a];
	  l[0][1] += vgm[0][1] * VV[i][a];
	  l[0][2] += vgm[0][2] * VV[1][a];
	  
	  l[1][0] += vgm[1][0] * VV[i][a];
	  l[1][1] += vgm[1][1] * VV[i][a];
	  l[1][2] += vgm[1][2] * VV[i][a];
	  
	  l[2][0] += vgm[2][0] * VV[i][a];
	  l[2][1] += vgm[2][1] * VV[i][a];
	  l[2][2] += vgm[2][2] * VV[i][a];
	  
	}
      }
    }
  }else{		
    /* cartesian */
    for(k = 1; k <= ppts; k++){
      for(a = 1; a <= ends; a++){
	/* velocity gradient matrix is transpose of grad v
	
	     | d_x(vx) d_y(vx) d_z(vx) |
	     | d_x(vy) d_y(vy) d_z(vy) |
	     | d_x(vz) d_y(vz) d_z(vz) |
	*/
	l[0][0] += GNx->ppt[GNPXINDEX(0, a, k)] * VV[1][a]; /* other contributions are zero */
	l[0][1] += GNx->ppt[GNPXINDEX(1, a, k)] * VV[1][a];
	l[0][2] += GNx->ppt[GNPXINDEX(2, a, k)] * VV[1][a];

	l[1][0] += GNx->ppt[GNPXINDEX(0, a, k)] * VV[2][a];
	l[1][1] += GNx->ppt[GNPXINDEX(1, a, k)] * VV[2][a];
	l[1][2] += GNx->ppt[GNPXINDEX(2, a, k)] * VV[2][a];

	l[2][0] += GNx->ppt[GNPXINDEX(0, a, k)] * VV[3][a];
	l[2][1] += GNx->ppt[GNPXINDEX(1, a, k)] * VV[3][a];
	l[2][2] += GNx->ppt[GNPXINDEX(2, a, k)] * VV[3][a];

      }
    }
  }
  if(ppts != 1){
    for(i=0;i<3;i++)
      for(j=0;j<3;j++)
	l[i][j] /= (float)ppts;
  }

}

void get_element_temperature(struct All_variables *E, float *ETEMP, int m)
{
  int e,kk;
  float TT[9];

  const int ends = enodes[E->mesh.nsd];
  const int nel = E->lmesh.nel;

  for(e = 1; e <= nel; e++)
  {
     ETEMP[e] = 0.0;
     for(kk = 1; kk <= ends; kk++)
     {
        TT[kk] = E->T[m][E->ien[m][e].node[kk]]; 
        ETEMP[e] += TT[kk];
     }
     ETEMP[e] /= ends;
  }
}

void get_element_yield(struct All_variables *E, float *EYIELD, int m)
{
  int e,kk;
  float YIE[9];

  const int ends = enodes[E->mesh.nsd];
  const int nel = E->lmesh.nel;

  for(e = 1; e <= nel; e++)
  {
     EYIELD[e] = 0.0;
     for(kk = 1; kk <= ends; kk++)
     {
        YIE[kk] = E->VIY[E->mesh.levmax][m][E->ien[m][e].node[kk]];
        EYIELD[e] += YIE[kk];
     }
     EYIELD[e] /= ends;
  }
}