[cig-commits] r17900 - in mc/1D/hc/trunk: . anisotropic_viscosity
becker at geodynamics.org
becker at geodynamics.org
Fri Feb 18 09:40:48 PST 2011
Author: becker
Date: 2011-02-18 09:40:47 -0800 (Fri, 18 Feb 2011)
New Revision: 17900
Added:
mc/1D/hc/trunk/anisotropic_viscosity/
mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c
mc/1D/hc/trunk/anisotropic_viscosity/anisotropic_viscosity.h
Modified:
mc/1D/hc/trunk/ggrd_grdtrack_util.c
Log:
Improved debugging output, and added anisotropic_viscosity
This will be used by both CitcomCU and CitcomS
Added: mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c
===================================================================
--- mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c (rev 0)
+++ mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c 2011-02-18 17:40:47 UTC (rev 17900)
@@ -0,0 +1,1298 @@
+/*
+
+set of routines to deal with anisotropic viscosity
+
+orthotropic viscosity following Muehlhaus, Moresi, Hobbs and Dufour
+(PAGEOPH, 159, 2311, 2002)
+
+tranverse isotropy following Han and Wahr (PEPI, 102, 33, 1997)
+
+
+*/
+
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+
+#include <math.h>
+#include "element_definitions.h"
+#include "global_defs.h"
+
+#ifdef CitcomS_global_defs_h /* CitcomS */
+
+#include "material_properties.h"
+#ifdef USE_GGRD
+#include "ggrd_handling.h"
+#endif
+
+
+#else /* CU */
+ void ggrd_read_anivisc_from_file(struct All_variables *);
+#endif
+
+#include "anisotropic_viscosity.h"
+void calc_cbase_at_tp(float , float , float *);
+void calc_cbase_at_tp_d(double , double , double *);
+#define CITCOM_DELTA(i,j) ((i==j)?(1.0):(0.0))
+
+void myerror_s(char *,struct All_variables *); /* for compatibility with CitcomS */
+
+/*
+
+output: D[6][6]
+
+input: n[3] director
+vis2: anisotropy factor
+avmode: anisotropy mode
+convert_to_spherical: 1/0 depending on geometry
+theta, phi : coordinates for conversion
+
+
+*/
+void get_constitutive(double D[6][6],double theta, double phi,
+ int convert_to_spherical,
+ float nx,float ny,float nz,
+ float vis2,
+ int avmode,
+ struct All_variables *E)
+{
+ double n[3];
+ if(E->viscosity.allow_anisotropic_viscosity){
+ if((E->monitor.solution_cycles == 0)&&
+ (E->viscosity.anivisc_start_from_iso)&&(E->monitor.visc_iter_count == 0)){
+ /*
+ first iteration of "stress-dependent" loop for first timestep
+ */
+ get_constitutive_isotropic(D);
+ }else{
+ /*
+ allow for a possibly anisotropic viscosity
+ */
+ n[0] = nx;n[1] = ny;n[2] = nz; /* Cartesian directors */
+ switch(avmode){
+ case CITCOM_ANIVISC_ORTHO_MODE:
+ /*
+ orthotropic
+ */
+ get_constitutive_orthotropic_viscosity(D,vis2,n,convert_to_spherical,theta,phi);
+ break;
+ case CITCOM_ANIVISC_TI_MODE:
+ /*
+ transversely isotropic
+ */
+ get_constitutive_ti_viscosity(D,vis2,0.,n,convert_to_spherical,theta,phi);
+ break;
+ default:
+ fprintf(stderr,"avmode %i\n",avmode);
+ myerror_s("get_constitutive: error, avmode undefined",E);
+ break;
+ }
+ //if(vis2 != 0) print_6x6_mat(stderr,D);
+ }
+ }else{
+ get_constitutive_isotropic(D);
+ }
+}
+
+
+/*
+
+transversely isotropic viscosity following Han and Wahr
+
+
+\nu_1 = isotropic viscosity, applies for e_31, e_23
+\nu_2 = weak anisotropy, applies for e_31, e_32
+\eta_1 = normal viscosity, (\eta_1+2\nu_1) control e_11, e_22
+\eta_2 = normal viscosity, (\eta_2+2\nu_2) = 2\eta_1 + 2\nu_1, controls e_33
+
+we use (for consistency with anisotropic viscosity)
+
+Delta = 1-\nu_2/\nu_1
+
+and
+
+\Gamma, such that \eta_1 = \Gamma \nu_1
+
+\nu_1 is the reference, isotropic viscosity, set to unity here, i.e.
+
+\nu_2 = 1 - \Delta ; \eta_1 = \Gamma ; (\eta_2 = 2 (\Gamma-\Delta)); for isotropy \Delta = 0, \Gamma = 0
+
+n[3] is the cartesian direction into which the weak shear points
+(ie. routine will rotate the 3 axis into the n direction) and will
+normalize n, if not already normalized
+
+
+*/
+ void get_constitutive_ti_viscosity(double D[6][6], double delta_vis, double gamma_vis,
+ double n[3], int convert_to_spherical,
+ double theta, double phi)
+{
+ double nlen,delta_vis2;
+ int ani;
+ /* isotropic part, in units of iso_visc */
+ get_constitutive_isotropic(D);
+ ani = FALSE;
+ if((fabs(delta_vis) > 3e-15) || (fabs(gamma_vis) > 3e-15)){
+ ani = TRUE;
+ /* get Cartesian anisotropy matrix by adding anisotropic
+ components */
+ D[0][0] += gamma_vis;
+ D[1][0] = D[0][1] = gamma_vis;
+ D[1][1] = D[0][0];
+ D[2][2] += 2.*gamma_vis;
+ D[4][4] -= delta_vis;
+ D[5][5] = D[4][4];
+ /*
+ the rotation routine takes care of normalization and will normalize n
+ */
+ //print_6x6_mat(stderr,D);
+ rotate_ti6x6_to_director(D,n); /* rotate such that the generic z
+ preferred axis is aligned with
+ the director */
+ //print_6x6_mat(stderr,D);fprintf(stderr,"\n");
+ }
+ if(ani && convert_to_spherical){
+ conv_cart6x6_to_spherical(D,theta,phi,D); /* rotate, can use in place */
+ }
+}
+
+ /*
+
+
+ compute a cartesian orthotropic anisotropic viscosity matrix (and
+ rotate it into CitcomS spherical, if requested)
+
+ viscosity is characterized by a eta_S (weak) viscosity in a shear
+ plane, to which the director is normal
+
+
+ output: D[0,...,5][0,...,5] constitutive matrix
+
+ input: delta_vis difference in viscosity from isotropic viscosity (set to unity here)
+
+ n[0,..,2]: director orientation, specify in cartesian
+
+
+ where delta_vis = (1 - eta_S/eta)
+
+
+ */
+ void get_constitutive_orthotropic_viscosity(double D[6][6], double delta_vis,
+ double n[3], int convert_to_spherical,
+ double theta, double phi)
+{
+ double nlen,delta_vis2;
+ double delta[3][3][3][3];
+ int ani;
+ ani=FALSE;
+ /* start with isotropic */
+ get_constitutive_isotropic(D);
+ /* get Cartesian anisotropy matrix */
+ if(fabs(delta_vis) > 3e-15){
+ ani = TRUE;
+ get_orth_delta(delta,n); /*
+ get anisotropy tensor, \Delta of
+ Muehlhaus et al. (2002)
+
+ this routine will normalize n, just in case
+
+ */
+ delta_vis2 = 2.0*delta_vis;
+ /* s_xx = s_tt */
+ D[0][0] -= delta_vis2 * delta[0][0][0][0]; /* * e_xx */
+ D[0][1] -= delta_vis2 * delta[0][0][1][1];
+ D[0][2] -= delta_vis2 * delta[0][0][2][2];
+ D[0][3] -= delta_vis * (delta[0][0][0][1]+delta[0][0][1][0]);
+ D[0][4] -= delta_vis * (delta[0][0][0][2]+delta[0][0][2][0]);
+ D[0][5] -= delta_vis * (delta[0][0][1][2]+delta[0][0][2][1]);
+
+ D[1][0] -= delta_vis2 * delta[1][1][0][0]; /* s_yy = s_pp */
+ D[1][1] -= delta_vis2 * delta[1][1][1][1];
+ D[1][2] -= delta_vis2 * delta[1][1][2][2];
+ D[1][3] -= delta_vis * (delta[1][1][0][1]+delta[1][1][1][0]);
+ D[1][4] -= delta_vis * (delta[1][1][0][2]+delta[1][1][2][0]);
+ D[1][5] -= delta_vis * (delta[1][1][1][2]+delta[1][1][2][1]);
+
+ D[2][0] -= delta_vis2 * delta[2][2][0][0]; /* s_zz = s_rr */
+ D[2][1] -= delta_vis2 * delta[2][2][1][1];
+ D[2][2] -= delta_vis2 * delta[2][2][2][2];
+ D[2][3] -= delta_vis * (delta[2][2][0][1]+delta[2][2][1][0]);
+ D[2][4] -= delta_vis * (delta[2][2][0][2]+delta[2][2][2][0]);
+ D[2][5] -= delta_vis * (delta[2][2][1][2]+delta[2][2][2][1]);
+
+ D[3][0] -= delta_vis2 * delta[0][1][0][0]; /* s_xy = s_tp */
+ D[3][1] -= delta_vis2 * delta[0][1][1][1];
+ D[3][2] -= delta_vis2 * delta[0][1][2][2];
+ D[3][3] -= delta_vis * (delta[0][1][0][1]+delta[0][1][1][0]);
+ D[3][4] -= delta_vis * (delta[0][1][0][2]+delta[0][1][2][0]);
+ D[3][5] -= delta_vis * (delta[0][1][1][2]+delta[0][1][2][1]);
+
+ D[4][0] -= delta_vis2 * delta[0][2][0][0]; /* s_xz = s_tr */
+ D[4][1] -= delta_vis2 * delta[0][2][1][1];
+ D[4][2] -= delta_vis2 * delta[0][2][2][2];
+ D[4][3] -= delta_vis * (delta[0][2][0][1]+delta[0][2][1][0]);
+ D[4][4] -= delta_vis * (delta[0][2][0][2]+delta[0][2][2][0]);
+ D[4][5] -= delta_vis * (delta[0][2][1][2]+delta[0][2][2][1]);
+
+ D[5][0] -= delta_vis2 * delta[1][2][0][0]; /* s_yz = s_pr */
+ D[5][1] -= delta_vis2 * delta[1][2][1][1];
+ D[5][2] -= delta_vis2 * delta[1][2][2][2];
+ D[5][3] -= delta_vis * (delta[1][2][0][1]+delta[1][2][1][0]);
+ D[5][4] -= delta_vis * (delta[1][2][0][2]+delta[1][2][2][0]);
+ D[5][5] -= delta_vis * (delta[1][2][1][2]+delta[1][2][2][1]);
+ }
+ if(ani && convert_to_spherical){
+ //print_6x6_mat(stderr,D);
+ conv_cart6x6_to_spherical(D,theta,phi,D); /* rotate, can use same mat for 6x6 */
+ //print_6x6_mat(stderr,D);fprintf(stderr,"\n");
+ }
+}
+void get_constitutive_isotropic(double D[6][6])
+{
+ /* isotropic part, in units of iso_visc */
+ zero_6x6(D);
+ D[0][0] = 2.0; /* xx = tt*/
+ D[1][1] = 2.0; /* yy = pp */
+ D[2][2] = 2.0; /* zz = rr */
+ D[3][3] = 1.0; /* xy = tp */
+ D[4][4] = 1.0; /* xz = rt */
+ D[5][5] = 1.0; /* yz = rp */
+}
+
+
+void set_anisotropic_viscosity_at_element_level(struct All_variables *E, int init_stage)
+{
+ int i,j,k,l,m,off,p,nel,elx,ely,elz,inode,elxlz,el,ani_layer;
+ double vis2,n[3],u,v,s,r,xloc[3],z_top,z_bottom;
+ float base[9],rout[3];
+ const int vpts = vpoints[E->mesh.nsd];
+ const int ends = enodes[E->mesh.nsd];
+ int mgmin,mgmax;
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ mgmin = E->mesh.gridmin;
+ mgmax = E->mesh.gridmax;
+#else /* Citcom CU */
+ mgmin = E->mesh.levmin;
+ mgmax = E->mesh.levmax;
+#endif
+ if(init_stage){
+ if(E->parallel.me == 0)
+ fprintf(stderr,"set_anisotropic_viscosity: allowing for %s viscosity\n",
+ (E->viscosity.allow_anisotropic_viscosity == 1)?("orthotropic"):("transversely isotropic"));
+ if(E->viscosity.anisotropic_viscosity_init)
+ myerror_s("anisotropic viscosity should not be initialized twice?!",E);
+ /* first call */
+ /* initialize anisotropic viscosity at element level, nodes will
+ get assigned later */
+ switch(E->viscosity.anisotropic_init){
+ case 0: /* isotropic */
+ case 3: /* first init for vel */
+ case 4: /* ISA */
+ case 5: /* same for mixed alignment */
+ if(E->parallel.me == 0)fprintf(stderr,"set_anisotropic_viscosity_at_element_level: initializing isotropic viscosity\n");
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ for(i=mgmin;i <= mgmax;i++){
+ nel = E->lmesh.NEL[i];
+ for (m=1;m <= E->sphere.caps_per_proc;m++) {
+ for(k=1;k <= nel;k++){
+ for(l=1;l <= vpts;l++){ /* assign to all integration points */
+ off = (k-1)*vpts + l;
+ E->EVI2[i][m][off] = 0.0;
+ E->EVIn1[i][m][off] = 1.0; E->EVIn2[i][m][off] = E->EVIn3[i][m][off] = 0.0;
+ E->avmode[i][m][off] = (unsigned char)
+ E->viscosity.allow_anisotropic_viscosity;
+ }
+ }
+ }
+ }
+#else /* CitcomCU */
+ for(i=mgmin;i <= mgmax;i++){
+ nel = E->lmesh.NEL[i];
+ for(k=1;k <= nel;k++){
+ for(l=1;l <= vpts;l++){ /* assign to all integration points */
+ off = (k-1)*vpts + l;
+ E->EVI2[i][off] = 0.0;
+ E->EVIn1[i][off] = 1.0; E->EVIn2[i][off] = E->EVIn3[i][off] = 0.0;
+ E->avmode[i][off] = (unsigned char)
+ E->viscosity.allow_anisotropic_viscosity;
+ }
+ }
+ }
+#endif
+ /* isotropic init end */
+ break;
+ case 1:
+ /*
+ random fluctuations, for testing a
+ worst case scenario
+
+ */
+ if(E->parallel.me == 0)fprintf(stderr,"set_anisotropic_viscosity_at_element_level: initializing random viscosity\n");
+ for(i=mgmin;i <= mgmax;i++){
+ nel = E->lmesh.NEL[i];
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ for (m=1;m <= E->sphere.caps_per_proc;m++) {
+#endif
+ for(k=1;k <= nel;k++){
+ vis2 = 0.01+drand48()*0.99; /* random fluctuation */
+ /* get random vector */
+ do{
+ u = -1 + drand48()*2;v = -1 + drand48()*2;
+ s = u*u + v*v;
+ }while(s > 1);
+ r = 2.0 * sqrt(1.0-s );
+ n[0] = u * r; /* x */
+ n[1] = v * r; /* y */
+ n[2] = 2.0*s -1 ; /* z */
+ for(l=1;l <= vpts;l++){ /* assign to all integration points */
+ off = (k-1)*vpts + l;
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ E->EVI2[i][m][off] = vis2;E->EVIn1[i][m][off] = n[0]; E->EVIn2[i][m][off] = n[1];E->EVIn3[i][m][off] = n[2];
+ E->avmode[i][m][off] = (unsigned char)E->viscosity.allow_anisotropic_viscosity;
+#else
+ E->EVI2[i][off] = vis2; E->EVIn1[i][off] = n[0]; E->EVIn2[i][off] = n[1]; E->EVIn3[i][off] = n[2];
+ E->avmode[i][off] = (unsigned char)E->viscosity.allow_anisotropic_viscosity;
+#endif
+ }
+ }
+#ifdef CitcomS_global_defs_h /* CitcomS m loop*/
+ }
+#endif
+ } /* mg loop */
+ /* end random init */
+ break;
+ case 2: /* from file */
+#ifndef USE_GGRD
+ fprintf(stderr,"set_anisotropic_viscosity_at_element_level: anisotropic_init mode 2 requires USE_GGRD compilation\n");
+ parallel_process_termination();
+#endif
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ if(E->sphere.caps == 12) /* global */
+ ggrd_read_anivisc_from_file(E,1);
+ else /* regional */
+ ggrd_read_anivisc_from_file(E,0);
+#else /* CitcomCU */
+ ggrd_read_anivisc_from_file(E);
+#endif
+ break;
+ case 6: /* tapered within layer */
+ if(E->parallel.me == 0)
+ fprintf(stderr,"set_anisotropic_viscosity_at_element_level: setting orthotropic tapered, vis2 min %g\n",
+ E->viscosity.ani_vis2_factor);
+ if(E->viscosity.anivisc_layer >= 0)myerror_s("set_anisotropic_viscosity_at_element_level: need to select layer",E);
+ ani_layer = -E->viscosity.anivisc_layer;
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ z_bottom = E->viscosity.zbase_layer[ani_layer-1];
+ if(ani_layer == 1)
+ z_top = E->sphere.ro;
+ else
+ z_top = E->viscosity.zbase_layer[ani_layer-2];
+#else /* CU */
+ z_bottom = E->viscosity.zbase_layer[ani_layer-1];
+ if(ani_layer == 1)
+ z_top = E->segment.zzlayer[E->segment.zlayers-1];
+ else
+ z_top = E->viscosity.zbase_layer[ani_layer-2];
+#endif
+ for(i=mgmin;i <= mgmax;i++){
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ for(m=1;m <= E->sphere.caps_per_proc;m++){
+#endif
+ elx = E->lmesh.ELX[i];elz = E->lmesh.ELZ[i];ely = E->lmesh.ELY[i];
+ elxlz = elx * elz;
+ for (j=1;j <= elz;j++){
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ if(E->mat[m][j] == -E->viscosity.anivisc_layer){
+#else
+ if(E->mat[j] == -E->viscosity.anivisc_layer){
+#endif
+ for(u=0.,inode=1;inode <= ends;inode++){ /* mean vertical coordinate */
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ off = E->ien[m][j].node[inode];
+ u += E->sx[m][3][off];
+#else
+ off = E->ien[j].node[inode];
+ if(E->control.Rsphere)
+ u += E->SX[3][off];
+ else
+ u += E->X[3][off];
+#endif
+ }
+ u /= ends;
+ /* do a log scale decrease of vis2 to ani_vis2_factor from bottom to top of layer */
+ vis2 = exp(log(E->viscosity.ani_vis2_factor) * (u-z_bottom)/(z_top-z_bottom));
+ //fprintf(stderr,"z %g (%g/%g) vis2 %g vis2_o %g frac %g\n",u,z_top,z_bottom,vis2, E->viscosity.ani_vis2_factor,(u-z_bottom)/(z_top-z_bottom));
+ /* 1-eta_s/eta */
+ vis2 = 1 - vis2;
+ for (k=1;k <= ely;k++){
+ for (l=1;l <= elx;l++) {
+ /* eq.(1) */
+ el = j + (l-1) * elz + (k-1)*elxlz;
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ xloc[0] = xloc[1] = xloc[2] = 0.0;
+ for(inode=1;inode <= ends;inode++){
+ off = E->ien[m][el].node[inode];
+ rtp2xyz((float)E->sx[m][3][off],(float)E->sx[m][1][off],(float)E->sx[m][2][off],rout);
+ xloc[0] += rout[0];xloc[1] += rout[1];xloc[2] += rout[2];
+ }
+ xloc[0]/=ends;xloc[1]/=ends;xloc[2]/=ends;xyz2rtp(xloc[0],xloc[1],xloc[2],rout);
+ calc_cbase_at_tp(rout[1],rout[2],base);convert_pvec_to_cvec(1.,0.,0.,base,rout);
+ n[0]=rout[0];n[1]=rout[1];n[2]=rout[2];
+ for(p=1;p <= vpts;p++){ /* assign to all integration points */
+ off = (el-1)*vpts + p;
+ E->EVI2[i][m][off] = vis2;
+ E->EVIn1[i][m][off] = n[0]; E->EVIn2[i][m][off] = n[1];E->EVIn3[i][m][off] = n[2];
+ E->avmode[i][m][off] = CITCOM_ANIVISC_ORTHO_MODE;
+ }
+#else /* CU */
+ if(E->control.Rsphere){ /* director in r direction */
+ xloc[0] = xloc[1] = xloc[2] = 0.0;
+ for(inode=1;inode <= ends;inode++){
+ off = E->ien[el].node[inode];
+ rtp2xyz((float)E->SX[3][off],(float)E->SX[1][off],(float)E->SX[2][off],rout);
+ xloc[0] += rout[0];xloc[1] += rout[1];xloc[2] += rout[2];
+ }
+ xloc[0]/=ends;xloc[1]/=ends;xloc[2]/=ends;xyz2rtp(xloc[0],xloc[1],xloc[2],rout);
+ calc_cbase_at_tp(rout[1],rout[2],base);convert_pvec_to_cvec(1.,0.,0.,base,rout);
+ n[0]=rout[0];n[1]=rout[1];n[2]=rout[2];
+ }else{ /* director in z direction */
+ n[0] = 0.;n[1] = 0.;n[2] = 1.;
+ }
+ for(p=1;p <= vpts;p++){ /* assign to all integration points */
+ off = (el-1)*vpts + p;
+ E->EVI2[i][off] = vis2;
+ E->EVIn1[i][off] = n[0]; E->EVIn2[i][off] = n[1];E->EVIn3[i][off] = n[2];
+ E->avmode[i][off] = CITCOM_ANIVISC_ORTHO_MODE;
+ }
+#endif
+ }
+ }
+ }else{ /* outside layer = isotropic */
+ for (k=1;k <= ely;k++){
+ for (l=1;l <= elx;l++){
+ el = j + (l-1) * elz + (k-1)*elxlz;
+ for(p=1;p <= vpts;p++){ /* assign to all integration points */
+ off = (el-1)*vpts + p;
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ E->EVI2[i][m][off] = 0;E->EVIn1[i][m][off] = 1; E->EVIn2[i][m][off] = 0;E->EVIn3[i][m][off] = 0;E->avmode[i][m][off] = CITCOM_ANIVISC_ORTHO_MODE;
+#else
+ E->EVI2[i][off] = 0;E->EVIn1[i][off] = 1; E->EVIn2[i][off] = 0;E->EVIn3[i][off] = 0;E->avmode[i][off] = CITCOM_ANIVISC_ORTHO_MODE;
+#endif
+ }
+ }
+ }
+ }
+ }
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ } /* end m loop */
+#endif
+ } /* end multigrid */
+ break; /* */
+ default:
+ fprintf(stderr,"set_anisotropic_viscosity_at_element_level: anisotropic_init %i undefined\n",
+ E->viscosity.anisotropic_init);
+ parallel_process_termination();
+ break;
+ }
+ E->viscosity.anisotropic_viscosity_init = TRUE;
+ /* end initialization stage */
+ }else{
+ //if(E->parallel.me == 0)fprintf(stderr,"reassigning anisotropic viscosity, mode %i\n",E->viscosity.anisotropic_init);
+ if((E->monitor.solution_cycles > 0) || (E->monitor.visc_iter_count > 0)){
+ /* standard operation every time step */
+ switch(E->viscosity.anisotropic_init){
+ /* if flow has been computed, use director aligned with ISA */
+ case 3: /* vel lignment */
+ /* we have a velocity solution already */
+
+ align_director_with_ISA_for_element(E,CITCOM_ANIVISC_ALIGN_WITH_VEL);
+ break;
+ case 4: /* vel alignment */
+ if((E->monitor.solution_cycles > 0) || (E->monitor.visc_iter_count > 0))
+ align_director_with_ISA_for_element(E,CITCOM_ANIVISC_ALIGN_WITH_ISA);
+ break;
+ case 5: /* mixed alignment */
+ if((E->monitor.solution_cycles > 0) || (E->monitor.visc_iter_count > 0))
+ align_director_with_ISA_for_element(E,CITCOM_ANIVISC_MIXED_ALIGN);
+ break;
+ default: /* default, no further modification of
+ anisotropy */
+ break;
+ }
+ }
+ }
+}
+
+#ifdef CitcomS_global_defs_h /* CitcomS */
+void normalize_director_at_nodes(struct All_variables *E,float **n1,float **n2, float **n3, int lev)
+{
+ int n,m;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(n=1;n<=E->lmesh.NNO[lev];n++){
+ normalize_vec3(&(n1[m][n]),&(n2[m][n]),&(n3[m][n]));
+ }
+}
+void normalize_director_at_gint(struct All_variables *E,float **n1,float **n2, float **n3, int lev)
+{
+ int m,e,i,enode;
+ const int nsd=E->mesh.nsd;
+ const int vpts=vpoints[nsd];
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++) {
+ enode = (e-1)*vpts+i;
+ normalize_vec3(&(n1[m][enode]),&(n2[m][enode]),&( n3[m][enode]));
+ }
+}
+#else /* CU */
+void normalize_director_at_nodes(struct All_variables *E,float *n1,float *n2, float *n3, int lev)
+{
+ int n;
+ for(n=1;n<=E->lmesh.NNO[lev];n++){
+ normalize_vec3(&(n1[n]),&(n2[n]),&(n3[n]));
+ }
+}
+void normalize_director_at_gint(struct All_variables *E,float *n1,float *n2, float *n3, int lev)
+{
+ int e,i,off;
+ const int nsd=E->mesh.nsd;
+ const int vpts=vpoints[nsd];
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++) {
+ off = (e-1)*vpts+i;
+ normalize_vec3(&(n1[off]),&(n2[off]),&(n3[off]));
+ }
+}
+#endif
+/*
+
+convert cartesian fourth order tensor (input c) to spherical, CitcomS
+format (output p)
+
+c and p cannot be the same matrix
+
+1: t 2: p 3: r
+
+(E only passed for debugging)
+
+*/
+
+ void conv_cart4x4_to_spherical(double c[3][3][3][3], double theta, double phi, double p[3][3][3][3])
+{
+ double rot[3][3];
+ get_citcom_spherical_rot(theta,phi,rot);
+ rot_4x4(c,rot,p);
+}
+
+/* convert [6][6] (input c) in cartesian to citcom spherical (output
+ p)
+
+ c and p can be the same amtrix
+
+*/
+void conv_cart6x6_to_spherical(double c[6][6],
+ double theta, double phi, double p[6][6])
+{
+ double c4[3][3][3][3],p4[3][3][3][3],rot[3][3];
+ get_citcom_spherical_rot(theta,phi,rot);
+ c4fromc6(c4,c);
+ rot_4x4(c4,rot,p4);
+ c6fromc4(p,p4);
+}
+/*
+
+rotate 6x6 D matrix with preferred axis aligned with z to the
+Cartesian director orientation, in place
+
+n will be normalized, just in case
+
+*/
+void rotate_ti6x6_to_director(double D[6][6],double n[3])
+{
+ double a[3][3][3][3],b[3][3][3][3],rot[3][3],
+ hlen2,x2,y2,xy,zm1;
+ /* normalize */
+ normalize_vec3d((n+0),(n+1),(n+2));
+ /* calc aux variable */
+ x2 = n[0]*n[0];y2 = n[1]*n[1];xy = n[0]*n[1];
+ hlen2 = x2 + y2;zm1 = n[2]-1;
+ if(hlen2 > 3e-15){
+ /* rotation matrix to get {0,0,1} to {x,y,z} */
+ rot[0][0] = (y2 + x2*n[2])/hlen2;
+ rot[0][1] = (xy*zm1)/hlen2;
+ rot[0][2] = n[0];
+ rot[1][0] = rot[0][1];
+ rot[1][1] = (x2 + y2*n[2])/hlen2;
+ rot[1][2] = n[1];
+ rot[2][0] = -n[0];
+ rot[2][1] = -n[1];
+ rot[2][2] = n[2];
+
+ /* rotate the D matrix */
+ c4fromc6(a,D);
+ rot_4x4(a,rot,b);
+ c6fromc4(D,b);
+ } /* already oriented right */
+
+}
+
+void get_citcom_spherical_rot(double theta, double phi, double rot[3][3]){
+ float base[9];
+ calc_cbase_at_tp((float)theta,(float)phi, base); /* compute cartesian basis at
+ theta, phi location */
+ rot[0][0] = base[3];rot[0][1] = base[4];rot[0][2] = base[5]; /* theta */
+ rot[1][0] = base[6];rot[1][1] = base[7];rot[1][2] = base[8]; /* phi */
+ rot[2][0] = base[0];rot[2][1] = base[1];rot[2][2] = base[2]; /* r */
+ //fprintf(stderr,"%g %g ; %g %g %g ; %g %g %g ; %g %g %g\n\n",
+ //theta,phi,rot[0][0],rot[0][1],rot[0][2],rot[1][0],rot[1][1],rot[1][2],rot[2][0],rot[2][1],rot[2][2]);
+}
+/*
+
+
+get fourth order anisotropy tensor for orthotropic viscosity from
+Muehlhaus et al. (2002)
+
+*/
+void get_orth_delta(double d[3][3][3][3],double n[3])
+{
+ int i,j,k,l;
+ double tmp;
+ normalize_vec3d((n+0),(n+1),(n+2));
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++)
+ for(k=0;k<3;k++)
+ for(l=0;l<3;l++){ /* eq. (4) from Muehlhaus et
+ al. (2002) */
+ tmp = n[i]*n[k]*CITCOM_DELTA(l,j);
+ tmp += n[j]*n[k]*CITCOM_DELTA(i,l);
+ tmp += n[i]*n[l]*CITCOM_DELTA(k,j);
+ tmp += n[j]*n[l]*CITCOM_DELTA(i,k);
+ tmp /= 2.0;
+ tmp -= 2*n[i]*n[j]*n[k]*n[l];
+ d[i][j][k][l] = tmp;
+ }
+
+}
+
+/*
+ mode =
+ CITCOM_ANIVISC_ALIGN_WITH_VEL: align with velocity
+ CITCOM_ANIVISC_ALIGN_WITH_ISA: align with ISA
+ CITCOM_ANIVISC_MIXED_ALIGN: mixed alignment
+
+
+*/
+void align_director_with_ISA_for_element(struct All_variables *E,
+ int mode)
+{
+ double rtf[4][9];
+ double VV[4][9],lgrad[3][3],isa[3],evel[3];
+ int e,i,off,m;
+ float base[9],n[3];
+ static struct CC Cc;
+ static struct CCX Ccx;
+ const int dims = E->mesh.nsd;
+ const int ppts = ppoints[dims];
+ const int vpts = vpoints[dims];
+ const int ends = enodes[dims];
+ const int lev = E->mesh.levmax;
+ const int nel = E->lmesh.nel;
+ unsigned char avmode;
+ double vis2 ;
+ /* anisotropy maximum strength */
+ vis2 = 1. - E->viscosity.ani_vis2_factor; /* 1-eta_w/eta_s */
+
+ if(E->parallel.me == 0){
+ switch(mode){
+ case CITCOM_ANIVISC_ALIGN_WITH_VEL:
+ fprintf(stderr,"align_director_with_ISA_for_element: aligning, max ani %g, align with vel\n",
+ vis2);
+ break;
+ case CITCOM_ANIVISC_ALIGN_WITH_ISA:
+ fprintf(stderr,"align_director_with_ISA_for_element: aligning, max ani %g, align with ISA\n",
+ vis2);
+ break;
+ case CITCOM_ANIVISC_MIXED_ALIGN:
+ fprintf(stderr,"align_director_with_ISA_for_element: aligning, max ani %g, align with uniaxial/vel\n",
+ vis2);
+ break;
+ default:
+ fprintf(stderr,"align_director_with_ISA_for_element: mode %i undefined\n",mode);
+ myerror_s("",E);
+ }
+ }
+ for(e=1; e <= nel; e++) {
+#ifdef CitcomS_global_defs_h /* CitcomS */
+ for(m=1;m <= E->sphere.caps_per_proc;m++) {
+ if(((E->viscosity.anivisc_layer > 0)&&
+ (E->mat[m][e] <= E->viscosity.anivisc_layer))||
+ ((E->viscosity.anivisc_layer < 0)&&
+ (E->mat[m][e] == -E->viscosity.anivisc_layer))){
+ get_rtf_at_ppts(E, m, lev, e, rtf); /* pressure points */
+ //if((e-1)%E->lmesh.elz==0)
+ construct_c3x3matrix_el(E,e,&E->element_Cc,&E->element_Ccx,lev,m,1);
+ for(i = 1; i <= ends; i++){ /* velocity at element nodes */
+ off = E->ien[m][e].node[i];
+ VV[1][i] = E->sphere.cap[m].V[1][off];
+ VV[2][i] = E->sphere.cap[m].V[2][off];
+ VV[3][i] = E->sphere.cap[m].V[3][off];
+ }
+ /* calculate velocity gradient matrix */
+ get_vgm_p(VV,&(E->N),&(E->GNX[lev][m][e]),&E->element_Cc,
+ &E->element_Ccx,rtf,E->mesh.nsd,ppts,ends,TRUE,lgrad,
+ evel);
+ /* calculate the ISA axis and determine the type of
+ anisotropy */
+ avmode = calc_isa_from_vgm(lgrad,evel,e,isa,E,mode);
+ /*
+ convert for spherical (Citcom system) to Cartesian
+ */
+ calc_cbase_at_tp(rtf[1][1],rtf[2][1],base);
+ convert_pvec_to_cvec(isa[2],isa[0],isa[1],base,n);
+ /* assign to director for all vpoints */
+ for(i=1;i <= vpts;i++){
+ off = (e-1)*vpts + i;
+ E->avmode[lev][m][off] = avmode;
+ E->EVI2[lev][m][off] = vis2;
+ E->EVIn1[lev][m][off] = n[0];
+ E->EVIn2[lev][m][off] = n[1];
+ E->EVIn3[lev][m][off] = n[2];
+ }
+ } /* in layer */
+ } /* caps */
+#else
+ if(((E->viscosity.anivisc_layer > 0)&&
+ (E->mat[e] <= E->viscosity.anivisc_layer))||
+ ((E->viscosity.anivisc_layer < 0)&&
+ (E->mat[e] == -E->viscosity.anivisc_layer))){
+ if(E->control.Rsphere){ /* need rtf for spherical */
+ get_rtf(E, e, 1, rtf, lev); /* pressure points */
+ //if((e-1)%E->lmesh.elz==0)
+ construct_c3x3matrix_el(E,e,&Cc,&Ccx,lev,1);
+ }
+ for(i = 1; i <= ends; i++){ /* velocity at element nodes */
+ off = E->ien[e].node[i];
+ VV[1][i] = E->V[1][off];
+ VV[2][i] = E->V[2][off];
+ VV[3][i] = E->V[3][off];
+ }
+ /* calculate velocity gradient matrix */
+ get_vgm_p(VV,&(E->N),&(E->GNX[lev][e]),&Cc, &Ccx,rtf,
+ E->mesh.nsd,ppts,ends,(E->control.Rsphere),lgrad,
+ evel);
+ /* calculate the ISA axis and determine the type of
+ anisotropy */
+ avmode = calc_isa_from_vgm(lgrad,evel,e,isa,E,mode);
+ /*
+ convert for spherical (Citcom system) to Cartesian
+ */
+ if(E->control.Rsphere){
+ calc_cbase_at_tp(rtf[1][1],rtf[2][1],base);
+ convert_pvec_to_cvec(isa[2],isa[0],isa[1],base,n);
+ }else{
+ n[0]=isa[0];n[1]=isa[1];n[2]=isa[2];
+ }
+ /* assign to director for all vpoints */
+ for(i=1;i <= vpts;i++){
+ off = (e-1)*vpts + i;
+ E->avmode[lev][off] = avmode;
+ E->EVI2[lev][off] = vis2;
+ E->EVIn1[lev][off] = n[0];
+ E->EVIn2[lev][off] = n[1];
+ E->EVIn3[lev][off] = n[2];
+ }
+ } /* end in layer branch */
+#endif
+ }
+
+}
+
+
+/*
+ compute the ISA axis from velocity gradient tensor l, element
+ velocity evel, and element number e
+
+ mode input:
+ CITCOM_ANIVISC_ALIGN_WITH_VEL: align with velocity
+ CITCOM_ANIVISC_ALIGN_WITH_ISA: align with ISA
+ CITCOM_ANIVISC_MIXED_ALIGN: mixed alignment
+
+ returns type of anisotropy CITCOM_ANIVISC_ORTHO_MODE : orthotropic (shear/normal)
+ CITCOM_ANIVISC_TI_MODE : TI
+
+*/
+unsigned char calc_isa_from_vgm(double l[3][3], double ev[3],
+ int e,double isa[3],
+ struct All_variables *E,
+ int mode)
+{
+ double d[3][3],r[3][3],ltrace,eval[3],lc[3][3],
+ evec[3][3],strain_scale,gol,t1,t2;
+ int i,j,isa_mode;
+ /* copy */
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++)
+ lc[i][j] = l[i][j];
+ remove_trace_3x3(lc);
+ calc_strain_from_vgm(lc,d); /* strain-rate */
+#ifndef USE_GGRD
+ myerror_s("need USE_GGRD compile for ISA axes",E);
+#else
+ ggrd_solve_eigen3x3(d,eval,evec,E); /* compute the eigensystem */
+#endif
+ /* normalize by largest abs(eigenvalue) */
+ strain_scale = ((t1=fabs(eval[2])) > (t2=fabs(eval[0])))?(t1):(t2);
+ for(i=0;i<3;i++){
+ eval[i] /= strain_scale; /* normalized eigenvalues */
+ for(j=0;j<3;j++)
+ lc[i][j] /= strain_scale;
+ }
+ /* recompute normalized strain rate and rotation */
+ calc_strain_from_vgm(lc,d);
+ calc_rot_from_vgm(lc,r); /* rotation */
+ switch(mode){
+ case CITCOM_ANIVISC_ALIGN_WITH_VEL: /* use velocity */
+ isa[0] = ev[0];isa[1] = ev[1];isa[2]=ev[2];
+ normalize_vec3d(isa,(isa+1),(isa+2));
+ return CITCOM_ANIVISC_TI_MODE;
+ break;
+ case CITCOM_ANIVISC_ALIGN_WITH_ISA:
+ isacalc(lc,&gol,isa,E,&isa_mode);
+ if((isa_mode == -1)||(isa_mode == 0)){
+ /* ISA cannot be found = align with flow */
+ isa[0] = ev[0];isa[1] = ev[1];isa[2]=ev[2];
+ normalize_vec3d(isa,(isa+1),(isa+2));
+ return CITCOM_ANIVISC_TI_MODE;
+ }else{ /* actual ISA */
+ return CITCOM_ANIVISC_ORTHO_MODE;
+ }
+ break;
+ case CITCOM_ANIVISC_MIXED_ALIGN:
+ /* mixed */
+ if(is_pure_shear(lc,d,r)){
+ /* align any pure shear state with the biggest absolute
+ eigenvector */
+ /* largest EV (extensional) */
+ //isa[0] = evec[0][0];isa[1] = evec[0][1];isa[2] = evec[0][2];
+ /* smallest (compressive) EV */
+ isa[0] = evec[2][0];isa[1] = evec[2][1];isa[2] = evec[2][2];
+ return CITCOM_ANIVISC_ORTHO_MODE;
+ }else{
+ /* simple shear */
+ isa[0] = ev[0];isa[1] = ev[1];isa[2]=ev[2]; /* align with vel for now */
+ normalize_vec3d(isa,(isa+1),(isa+2));
+ return CITCOM_ANIVISC_TI_MODE; /* TI */
+ }
+ break;
+ default:
+ myerror_s("ISA mode undefined",E);
+ break;
+ }
+ return 0;
+}
+/*
+ determine if deformation state is pure shear from normalized
+ velocity gradient, strain, and rotation tensor
+
+*/
+
+int is_pure_shear(double l[3][3],double e[3][3],double r[3][3])
+{
+
+ double mrot,tmp,e2;
+ /* find max rotation */
+ mrot = fabs(r[0][1]); /* xy */
+ if((tmp=fabs(r[0][2]))>mrot)mrot = tmp; /* xz */
+ if((tmp=fabs(r[1][2]))>mrot)mrot = tmp; /* yz */
+ /* second invariant of strain-rate tensor */
+ e2 = second_invariant_from_3x3(e);
+ if(mrot/e2 >= 1)
+ return 0; /* simple shear */
+ else
+ return 1; /* pure shear */
+
+}
+/*
+ rotate fourth order tensor
+ c4 and c4c cannot be the same matrix
+
+*/
+void rot_4x4(double c4[3][3][3][3], double r[3][3],
+ double c4c[3][3][3][3])
+{
+
+ int i1,i2,i3,i4,j1,j2,j3,j4;
+
+ zero_4x4(c4c);
+
+ for(i1=0;i1<3;i1++)
+ for(i2=0;i2<3;i2++)
+ for(i3=0;i3<3;i3++)
+ for(i4=0;i4<3;i4++)
+ for(j1=0;j1<3;j1++)
+ for(j2=0;j2<3;j2++)
+ for(j3=0;j3<3;j3++)
+ for(j4=0;j4<3;j4++)
+ c4c[i1][i2][i3][i4] +=
+ r[i1][j1] * r[i2][j2] *
+ r[i3][j3] * r[i4][j4] * c4[j1][j2][j3][j4];
+
+}
+void zero_6x6(double a[6][6])
+{
+ int i,j;
+ for(i=0;i<6;i++)
+ for(j=0;j<6;j++)
+ a[i][j] = 0.;
+
+}
+void zero_4x4(double a[3][3][3][3])
+{
+ int i1,i2,i3,i4;
+ for(i1=0;i1<3;i1++)
+ for(i2=0;i2<3;i2++)
+ for(i3=0;i3<3;i3++)
+ for(i4=0;i4<3;i4++)
+ a[i1][i2][i3][i4] = 0.0;
+
+}
+void copy_4x4(double a[3][3][3][3], double b[3][3][3][3])
+{
+
+ int i1,i2,i3,i4;
+ for(i1=0;i1<3;i1++)
+ for(i2=0;i2<3;i2++)
+ for(i3=0;i3<3;i3++)
+ for(i4=0;i4<3;i4++)
+ b[i1][i2][i3][i4] = a[i1][i2][i3][i4];
+}
+void copy_6x6(double a[6][6], double b[6][6])
+{
+
+ int i1,i2;
+ for(i1=0;i1<6;i1++)
+ for(i2=0;i2<6;i2++)
+ b[i1][i2] = a[i1][i2];
+}
+
+void print_6x6_mat(FILE *out, double c[6][6])
+{
+ int i,j;
+ for(i=0;i<6;i++){
+ for(j=0;j<6;j++)
+ fprintf(out,"%14.5e ",(fabs(c[i][j])<5e-15)?(0):(c[i][j]));
+ fprintf(out,"\n");
+ }
+}
+void print_3x3_mat(FILE *out, double c[3][3])
+{
+ int i,j;
+ for(i=0;i<3;i++){
+ for(j=0;j<3;j++)
+ fprintf(out,"%14.5e ",(fabs(c[i][j])<5e-15)?(0):(c[i][j]));
+ fprintf(out,"\n");
+ }
+}
+/*
+ create a fourth order tensor representation from the voigt
+ notation, assuming only upper half is filled in
+
+*/
+void c4fromc6(double c4[3][3][3][3],double c[6][6])
+{
+ int i,j;
+
+ c4[0][0][0][0] = c[0][0];
+ c4[0][0][1][1] = c[0][1];
+ c4[0][0][2][2] = c[0][2];
+ c4[0][0][0][1] = c4[0][0][1][0] = c[0][3];
+ c4[0][0][0][2] = c4[0][0][2][0] = c[0][4];
+ c4[0][0][1][2] = c4[0][0][2][1] = c[0][5];
+
+ c4[1][1][0][0] = c[0][1];
+ c4[1][1][1][1] = c[1][1];
+ c4[1][1][2][2] = c[1][2];
+ c4[1][1][0][1] = c4[1][1][1][0] = c[1][3];
+ c4[1][1][0][2] = c4[1][1][2][0] = c[1][4];
+ c4[1][1][1][2] = c4[1][1][2][1] = c[1][5];
+
+ c4[2][2][0][0] = c[0][2];
+ c4[2][2][1][1] = c[1][2];
+ c4[2][2][2][2] = c[2][2];
+ c4[2][2][0][1] = c4[2][2][1][0] = c[2][3];
+ c4[2][2][0][2] = c4[2][2][2][0] = c[2][4];
+ c4[2][2][1][2] = c4[2][2][2][1] = c[2][5];
+
+ c4[0][1][0][0] = c[0][3];
+ c4[0][1][1][1] = c[1][3];
+ c4[0][1][2][2] = c[2][3];
+ c4[0][1][0][1] = c4[0][1][1][0] = c[3][3];
+ c4[0][1][0][2] = c4[0][1][2][0] = c[3][4];
+ c4[0][1][1][2] = c4[0][1][2][1] = c[3][5];
+
+ c4[0][2][0][0] = c[0][4];
+ c4[0][2][1][1] = c[1][4];
+ c4[0][2][2][2] = c[2][4];
+ c4[0][2][0][1] = c4[0][2][1][0] = c[3][4];
+ c4[0][2][0][2] = c4[0][2][2][0] = c[4][4];
+ c4[0][2][1][2] = c4[0][2][2][1] = c[4][5];
+
+ c4[1][2][0][0] = c[0][5];
+ c4[1][2][1][1] = c[1][5];
+ c4[1][2][2][2] = c[2][5];
+ c4[1][2][0][1] = c4[1][2][1][0] = c[3][5];
+ c4[1][2][0][2] = c4[1][2][2][0] = c[4][5];
+ c4[1][2][1][2] = c4[1][2][2][1] = c[5][5];
+
+ /* assign the symmetric diagonal terms */
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++){
+ c4[1][0][i][j] = c4[0][1][i][j];
+ c4[2][0][i][j] = c4[0][2][i][j];
+ c4[2][1][i][j] = c4[1][2][i][j];
+ }
+
+}
+void c6fromc4(double c[6][6],double c4[3][3][3][3])
+{
+ int i,j;
+
+ c[0][0] = c4[0][0][0][0];
+ c[0][1] = c4[0][0][1][1];
+ c[0][2] = c4[0][0][2][2];
+ c[0][3] = c4[0][0][0][1];
+ c[0][4] = c4[0][0][0][2];
+ c[0][5] = c4[0][0][1][2];
+
+ c[1][1] = c4[1][1][1][1];
+ c[1][2] = c4[1][1][2][2];
+ c[1][3] = c4[1][1][0][1];
+ c[1][4] = c4[1][1][0][2];
+ c[1][5] = c4[1][1][1][2];
+
+ c[2][2] = c4[2][2][2][2];
+ c[2][3] = c4[2][2][0][1];
+ c[2][4] = c4[2][2][0][2];
+ c[2][5] = c4[2][2][1][2];
+
+ c[3][3] = c4[0][1][0][1];
+ c[3][4] = c4[0][1][0][2];
+ c[3][5] = c4[0][1][1][2];
+
+ c[4][4] = c4[0][2][0][2];
+ c[4][5] = c4[0][2][1][2];
+
+ c[5][5] = c4[1][2][1][2];
+ /* fill in the lower half */
+ for(i=0;i<6;i++)
+ for(j=i+1;j<6;j++)
+ c[j][i] = c[i][j];
+}
+
+void isacalc(double l[3][3], double *gol,double isa[3],
+ struct All_variables *E,int *isa_mode)
+{
+
+ //
+ // input: l: *normalized* velocity gradient tensor, in FORTRAN sorting
+ // output: isa(3): infite strain axis,
+ // gol: grain orientation lag
+ double ltrace,isa_diff;
+ double f[3][3],eval[3],evec[3][3],u[3][3],le[3][3];
+
+ // make sure l does not have a trace
+ remove_trace_3x3(l);
+
+#ifdef CITCOM_USE_EXPOKIT
+ calc_exp_matrixt(l,75,le,E); /* following Kaminski & Ribe (G-Cubed, 2001) */
+ f_times_ft(le,u);ggrd_solve_eigen3x3(u,eval,evec,E);
+ isa[0] = evec[0][0]; isa[1] = evec[0][1]; isa[2] = evec[0][2];
+ calc_exp_matrixt(l,80,le,E); f_times_ft(le,u);ggrd_solve_eigen3x3(u,eval,evec,E);
+ isa_diff = 1.-fabs(evec[0][0]*isa[0]+evec[0][1]*isa[1]+evec[0][2]*isa[2]);
+ if(isa_diff > 1e-4) /* ISA does not exist */
+ *isa_mode = -1;
+ else
+ *isa_mode = 1;
+ //fprintf(stderr,"A: %11g %11g %11g - %11g %i\n",isa[0],isa[1],isa[2],isa_diff,*isa_mode);
+
+#else
+ // Limit deformation gradient tensor for infinite time
+ // calculation of the ISE orientation using Sylvester's formula
+ drex_eigen(l,f,isa_mode);
+ if(*isa_mode == -1){
+ // isa is flow1
+ isa[0]=isa[1]=isa[2] = -1.;
+ }else if(*isa_mode == 0){
+ isa[0] = isa[1] = isa[2] = 0;
+ *gol = -1.;
+ }else{
+ // 2. formation of the left-stretch tensor U = FFt
+ f_times_ft(f,u);
+ // 3. eigen-values and eigen-vectors of U
+ ggrd_solve_eigen3x3(u,eval,evec,E);
+ // largest eigenvector
+ isa[0] = evec[0][0];isa[1] = evec[0][1];isa[2] = evec[0][2];
+ //fprintf(stderr,"B: %11g %11g %11g\n",evec[0][0],evec[0][1],evec[0][2]);
+ }
+#endif
+}
+/*
+ F^2 = F * F^T
+*/
+void f_times_ft(double f[3][3],double out[3][3])
+{
+ int i,j,k;
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++){
+ out[i][j] = 0.0;
+ for(k=0;k<3;k++)
+ out[i][j] += f[i][k] * f[j][k];
+ }
+}
+/* find eigenvalues of velocity gradient tensor (modified from DREX
+ code of Kaminski et al. 2004)
+
+*/
+void drex_eigen(double l[3][3],double f[3][3], int *mode)
+{
+ double a2,a3,q,q3,r2,r,theta,xx,lambda1,lambda2,lambda3,sq2;
+ const double four_pi = 4.0*M_PI,
+ two_pi = 2.0*M_PI,
+ one_third=1./3.;
+ /* looking for the eigen-values of L (using tr(l)=0) */
+ a2 = l[0][0] * l[1][1] + l[1][1] * l[2][2] + l[2][2]*l[0][0] -
+ l[0][1] * l[1][0] - l[1][2]*l[2][1] - l[2][0]*l[0][2];
+ a3 = l[0][0]*l[1][2]*l[2][1] + l[0][1]*l[1][0]*l[2][2] +
+ l[0][2]*l[1][1]*l[2][0] - l[0][0]*l[1][1]*l[2][2] -
+ l[0][1]*l[1][2]*l[2][0] - l[0][2]*l[1][0]*l[2][1];
+
+ q = -a2/3.;
+ r = a3/2.;
+ q3 = q*q*q;
+ r2 = r*r;
+
+ if(fabs(q) < 1e-9){
+ /* simple shear, isa=veloc */
+ *mode = -1;
+ }else if(q3-r2 >= 0){
+ sq2 = 2*sqrt(q);
+
+ theta = acos(pow(r/q,1.5));
+ lambda1 = -sq2*cos(theta/3);
+ lambda2 = -sq2*cos((theta+two_pi)/3.);
+ lambda3 = -sq2*cos((theta+four_pi)/3.);
+
+ if (fabs(lambda1-lambda2) < 1e-13)
+ lambda1 = lambda2;
+ if (fabs(lambda2-lambda3) < 1e-13)
+ lambda2 = lambda3;
+ if (fabs(lambda3-lambda1) < 1e-13)
+ lambda3 = lambda1;
+
+ if((lambda1 > lambda2) && (lambda1 > lambda3)) {
+ malmul_scaled_id(f,l,-lambda2,-lambda3);*mode=1;
+ }else if((lambda2 > lambda3 ) && (lambda2 > lambda1)){
+ malmul_scaled_id(f,l,-lambda3,-lambda1);*mode = 1;
+ }else if((lambda3 > lambda1 ) && (lambda3 > lambda2)){
+ malmul_scaled_id(f,l,-lambda1,-lambda2);*mode = 1;
+ }else if((lambda1 == lambda2 ) && (lambda3 > lambda1)) {
+ malmul_scaled_id(f,l,-lambda1,-lambda2);*mode = 1;
+ }else if((lambda2 == lambda3 ) && (lambda1 > lambda2)){
+ malmul_scaled_id(f,l,-lambda2,-lambda3);*mode = 1;
+ }else if((lambda3 == lambda1 ) && (lambda2 > lambda3)){
+ malmul_scaled_id(f,l,-lambda3,-lambda1);*mode = 1;
+ }else if((lambda1 == lambda2 ) && (lambda3 < lambda1)){
+ *mode =0;
+ }else if((lambda2 == lambda3 ) && (lambda1 < lambda2)){
+ *mode = 0;
+ }else if((lambda3 == lambda1 ) && (lambda2 < lambda3)){
+ *mode = 0;
+ }
+ }else{
+ xx = pow(sqrt(r2-q3)+fabs(r),one_third);
+ if(r < 0)
+ lambda1 = xx+q/xx;
+ else
+ lambda1 = -xx+q/xx;
+ lambda2 = -lambda1/2.;
+ lambda3 = -lambda1/2.;
+ if (lambda1 > 1e-9) {
+ malmul_scaled_id(f,l,-lambda2,-lambda3);
+ *mode = 2;
+ }else{
+ *mode = 0;
+ }
+
+ }
+}
+/* f = matmul(l-lambda2*Id,l-lambda3*Id); */
+void malmul_scaled_id(double f[3][3],double l[3][3],
+ double f1,double f2)
+{
+ double a[3][3],b[3][3];
+ int i,j;
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++){
+ a[i][j] = b[i][j] = l[i][j];
+ if(i==j){
+ a[i][j] += f1;
+ b[i][j] += f2;
+ }
+ }
+ matmul_3x3(a,b,f);
+
+
+}
+#ifdef CITCOM_USE_EXPOKIT
+/*
+
+calculate exp(A t) using DGPADM from EXPOKIT for a 3x3
+
+(needs LAPACK)
+
+*/
+void calc_exp_matrixt(double a[3][3],double t,double ae[3][3],
+ struct All_variables *E)
+{
+ int ideg=6;// degre of Pade approximation, six should be ok
+ int m=3,ldh=3;// a(ldh,m) dimensions
+ double *wsp;// work space
+ double af[9];
+ int ipiv[4],iexph,ns;// workspace, output pointer, nr of squareas
+ int iflag,lwsp;// exit code, size of workspace
+ int i,j,k;
+ /*
+ work space
+ */
+ lwsp = 2*(4*m*m+ideg+1);// size of workspace, oversized by factor two
+ wsp = (double *)calloc(lwsp,sizeof(double));
+ /* assign fortran style */
+ for(k=i=0;i<3;i++)
+ for(j=0;j<3;j++,k++)
+ af[k] = a[i][j];
+ //
+ // call to expokit routine
+ dgpadm_(&ideg,&m,&t,af,&ldh,wsp,&lwsp,ipiv,&iexph,&ns,&iflag);
+ if(iflag < 0)
+ myerror_s("calc_exp_matrixt: problem in dgpadm",E);
+ // assign to output
+ for(i=0,k=iexph-1;i<3;i++)
+ for(j=0;j<3;j++,k++)
+ ae[i][j] = wsp[k];
+
+ free(wsp);
+
+}
+#endif
+#ifdef CitcomS_global_defs_h /* CitcomS */
+void myerror_s(char *a,struct All_variables *E){
+ myerror(E,a);
+}
+#else /* CU */
+void myerror_s(char *a,struct All_variables *E){
+ myerror(a,E);
+}
+#endif
+
+
+
+#endif
+
Added: mc/1D/hc/trunk/anisotropic_viscosity/anisotropic_viscosity.h
===================================================================
--- mc/1D/hc/trunk/anisotropic_viscosity/anisotropic_viscosity.h (rev 0)
+++ mc/1D/hc/trunk/anisotropic_viscosity/anisotropic_viscosity.h 2011-02-18 17:40:47 UTC (rev 17900)
@@ -0,0 +1,63 @@
+#ifndef __CITCOM_READ_ANIVISC_HEADER__
+
+#define CITCOM_ANIVISC_ORTHO_MODE 1
+#define CITCOM_ANIVISC_TI_MODE 2
+
+#define CITCOM_ANIVISC_ALIGN_WITH_VEL 0
+#define CITCOM_ANIVISC_ALIGN_WITH_ISA 1
+#define CITCOM_ANIVISC_MIXED_ALIGN 2
+
+void get_constitutive(double [6][6], double, double, int, float,float,float,float,int,struct All_variables *);
+void get_constitutive_ti_viscosity(double [6][6], double, double, double [3], int, double, double);
+void get_constitutive_orthotropic_viscosity(double [6][6], double, double [3], int, double, double);
+void get_constitutive_isotropic(double [6][6]);
+void set_anisotropic_viscosity_at_element_level(struct All_variables *, int);
+
+
+
+void conv_cart4x4_to_spherical(double [3][3][3][3], double, double, double [3][3][3][3]);
+void conv_cart6x6_to_spherical(double [6][6], double, double, double [6][6]);
+void rotate_ti6x6_to_director(double [6][6], double [3]);
+void get_citcom_spherical_rot(double, double, double [3][3]);
+void get_orth_delta(double [3][3][3][3], double [3]);
+void align_director_with_ISA_for_element(struct All_variables *, int);
+unsigned char calc_isa_from_vgm(double [3][3], double [3], int, double [3], struct All_variables *, int);
+int is_pure_shear(double [3][3], double [3][3], double [3][3]);
+void rot_4x4(double [3][3][3][3], double [3][3], double [3][3][3][3]);
+void zero_6x6(double [6][6]);
+void zero_4x4(double [3][3][3][3]);
+void copy_4x4(double [3][3][3][3], double [3][3][3][3]);
+void copy_6x6(double [6][6], double [6][6]);
+void print_6x6_mat(FILE *, double [6][6]);
+void c4fromc6(double [3][3][3][3], double [6][6]);
+void c6fromc4(double [6][6], double [3][3][3][3]);
+void isacalc(double [3][3], double *, double [3], struct All_variables *, int *);
+void f_times_ft(double [3][3], double [3][3]);
+void drex_eigen(double [3][3], double [3][3], int *);
+void malmul_scaled_id(double [3][3], double [3][3], double, double);
+
+
+void print_3x3_mat(FILE *, double [3][3]);
+
+void calc_exp_matrixt(double [3][3],double ,double [3][3],
+ struct All_variables *);
+
+void dgpadm_(int *,int *,double *,double *,int *,double *,int *,
+ int *,int *,int *,int *);
+void get_vgm_p(double [4][9],struct Shape_function *,
+ struct Shape_function_dx *,
+ struct CC *, struct CCX *, double [4][9],
+ int ,int , int , int ,
+ double [3][3], double [3]);
+
+#ifdef CitcomS_global_defs_h /* CitcomS */
+void normalize_director_at_nodes(struct All_variables *, float **, float **, float **, int);
+void normalize_director_at_gint(struct All_variables *, float **, float **, float **, int);
+#else
+void normalize_director_at_nodes(struct All_variables *, float *, float *, float *, int);
+void normalize_director_at_gint(struct All_variables *, float *, float *, float *, int);
+
+#endif
+
+#define __CITCOM_READ_ANIVISC_HEADER__
+#endif
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2011-02-18 03:33:07 UTC (rev 17899)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2011-02-18 17:40:47 UTC (rev 17900)
@@ -683,7 +683,7 @@
*f = (float *) calloc((*mm) * (*nz) ,sizeof (float));
if(!(*f)){
- fprintf(stderr,"ggrd_grdtrack_init: f memory error\n");
+ fprintf(stderr,"ggrd_grdtrack_init: f memory error, mm: %i (%i by %i) by nz: %i \n",*mm,mx,my, *nz);
return 9;
}
if(verbose >= 2){
More information about the CIG-COMMITS
mailing list