[cig-commits] r17904 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Fri Feb 18 11:45:49 PST 2011
Author: becker
Date: 2011-02-18 11:45:49 -0800 (Fri, 18 Feb 2011)
New Revision: 17904
Added:
mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
Removed:
mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
Modified:
mc/3D/CitcomS/trunk/lib/
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
Log:
Modified storage of anisotropic viscosity files. Those are now
Property changes on: mc/3D/CitcomS/trunk/lib
___________________________________________________________________
Name: svn:externals
+ anisotropic_viscosity http://geodynamics.org/svn/cig/mc/1D/hc/trunk/anisotropic_viscosity
Deleted: mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2011-02-18 19:44:31 UTC (rev 17903)
+++ mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2011-02-18 19:45:49 UTC (rev 17904)
@@ -1,1298 +0,0 @@
-/*
-
-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/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c (rev 0)
+++ mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2011-02-18 19:45:49 UTC (rev 17904)
@@ -0,0 +1 @@
+link anisotropic_viscosity/Anisotropic_viscosity.c
\ No newline at end of file
Property changes on: mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
___________________________________________________________________
Name: svn:special
+ *
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2011-02-18 19:44:31 UTC (rev 17903)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2011-02-18 19:45:49 UTC (rev 17904)
@@ -286,6 +286,7 @@
/* write nodal coordinate to file, big endian */
for(j=1;j <= E->sphere.caps_per_proc;j++) {
for(i=1;i <= E->lmesh.nno;i++) {
+ /* cartesian coordinates */
x[0]=E->x[j][1][i];x[1]=E->x[j][2][i];x[2]=E->x[j][3][i];
if(be_write_float_to_file(x,3,fp1) != 3)
BE_WERROR;
Deleted: mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h 2011-02-18 19:44:31 UTC (rev 17903)
+++ mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h 2011-02-18 19:45:49 UTC (rev 17904)
@@ -1,63 +0,0 @@
-#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
Added: mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h (rev 0)
+++ mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h 2011-02-18 19:45:49 UTC (rev 17904)
@@ -0,0 +1 @@
+link anisotropic_viscosity/anisotropic_viscosity.h
\ No newline at end of file
Property changes on: mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
___________________________________________________________________
Name: svn:special
+ *
More information about the CIG-COMMITS
mailing list