[cig-commits] r17212 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Thu Sep 23 20:35:36 PDT 2010
Author: becker
Date: 2010-09-23 20:35:36 -0700 (Thu, 23 Sep 2010)
New Revision: 17212
Modified:
mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Instructions.c
Log:
Augmented G matrix with anisotropic computation, may or may not be a good idea.
Modified: mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-22 23:48:53 UTC (rev 17211)
+++ mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-24 03:35:36 UTC (rev 17212)
@@ -52,7 +52,7 @@
void get_constitutive(double D[6][6], int lev, int m,
- int off, double theta, double phi,
+ int enode, double theta, double phi,
int convert_to_spherical,
struct All_variables *E)
{
@@ -64,15 +64,15 @@
get_constitutive_isotropic(D);
else{
/* allow for a possibly anisotropic viscosity */
- n[0] = E->EVIn1[lev][m][off];
- n[1] = E->EVIn2[lev][m][off];
- n[2] = E->EVIn3[lev][m][off]; /* Cartesian directors */
+ n[0] = E->EVIn1[lev][m][enode];
+ n[1] = E->EVIn2[lev][m][enode];
+ n[2] = E->EVIn3[lev][m][enode]; /* Cartesian directors */
if(E->viscosity.allow_anisotropic_viscosity == 1){ /* orthotropic */
- get_constitutive_orthotropic_viscosity(D,E->EVI2[lev][m][off],
+ get_constitutive_orthotropic_viscosity(D,E->EVI2[lev][m][enode],
n,convert_to_spherical,theta,phi);
}else if(E->viscosity.allow_anisotropic_viscosity == 2){
/* transversely isotropic */
- get_constitutive_ti_viscosity(D,E->EVI2[lev][m][off],0.,
+ get_constitutive_ti_viscosity(D,E->EVI2[lev][m][enode],0.,
n,convert_to_spherical,theta,phi);
}
}
@@ -247,7 +247,7 @@
}
void set_anisotropic_viscosity_at_element_level(struct All_variables *E, int init_stage)
{
- int i,j,k,l,off,nel;
+ int i,j,k,l,enode,nel;
double vis2,n[3],u,v,s,r;
const int vpts = vpoints[E->mesh.nsd];
@@ -266,9 +266,9 @@
for (j=1;j<=E->sphere.caps_per_proc;j++) {
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][j][off] = 0.0;
- E->EVIn1[i][j][off] = 1.0; E->EVIn2[i][j][off] = E->EVIn3[i][j][off] = 0.0;
+ enode = (k-1)*vpts + l;
+ E->EVI2[i][j][enode] = 0.0;
+ E->EVIn1[i][j][enode] = 1.0; E->EVIn2[i][j][enode] = E->EVIn3[i][j][enode] = 0.0;
}
}
}
@@ -299,11 +299,11 @@
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;
- E->EVI2[i][j][off] = vis2;
- E->EVIn1[i][j][off] = n[0];
- E->EVIn2[i][j][off] = n[1];
- E->EVIn3[i][j][off] = n[2];
+ enode = (k-1)*vpts + l;
+ E->EVI2[i][j][enode] = vis2;
+ E->EVIn1[i][j][enode] = n[0];
+ E->EVIn2[i][j][enode] = n[1];
+ E->EVIn3[i][j][enode] = n[2];
}
}
}
@@ -343,14 +343,14 @@
}
void normalize_director_at_gint(struct All_variables *E,float **n1,float **n2, float **n3, int lev)
{
- int m,e,i,off;
+ 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++) {
- off = (e-1)*vpts+i;
- normalize_vec3(&(n1[m][off]),&(n2[m][off]),&( n3[m][off]));
+ enode = (e-1)*vpts+i;
+ normalize_vec3(&(n1[m][enode]),&(n2[m][enode]),&( n3[m][enode]));
}
}
/*
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2010-09-22 23:48:53 UTC (rev 17211)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2010-09-24 03:35:36 UTC (rev 17212)
@@ -319,11 +319,13 @@
off = (el-1)*vpts+k;
W[k]=g_point[k].weight[dims-1]*E->GDA[lev][m][el].vpt[k]*E->EVI[lev][m][off];
#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
- /* allow for a possibly anisotropic viscosity */
- get_constitutive(D[k],lev,m,off,rtf[1][k],rtf[2][k],TRUE,E);
+ if(E->viscosity.allow_anisotropic_viscosity){
+ /* allow for a possibly anisotropic viscosity */
+ get_constitutive(D[k],lev,m,off,rtf[1][k],rtf[2][k],TRUE,E);
+ }
#endif
}
-
+ /* */
get_ba(&(E->N), &(E->GNX[lev][m][el]), &E->element_Cc, &E->element_Ccx,
rtf, E->mesh.nsd, ba);
@@ -906,6 +908,7 @@
/*==============================================================
Function to supply the element g matrix for a given element e.
+ used for the divergence
============================================================== */
void get_elt_g(E,el,elt_del,lev,m)
@@ -916,13 +919,20 @@
{
void get_rtf_at_ppts();
- int p,a,i;
+ int p,a,i,j,k;
double ra,ct,si,x[4],rtf[4][9];
double temp;
-
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ double Dtmp[6][6],Duse[6][6],rtf2[4][9],weight;
+ const int vpts=VPOINTS3D;
+ const int modify_g = TRUE;
+ //const int modify_g = FALSE;
+ int off;
+ double ba[9][9][4][7];
+#endif
const int dims=E->mesh.nsd;
const int ends=enodes[dims];
-
+
/* Special case, 4/8 node bilinear cartesian square/cube element -> 1 pressure point */
if ((el-1)%E->lmesh.ELZ[lev]==0)
@@ -930,31 +940,69 @@
get_rtf_at_ppts(E, m, lev, el, rtf);
- temp=p_point[1].weight[dims-1] * E->GDA[lev][m][el].ppt[1];
+ temp = p_point[1].weight[dims-1] * E->GDA[lev][m][el].ppt[1];
- ra = rtf[3][1];
- si = 1.0/sin(rtf[1][1]);
- ct = cos(rtf[1][1])*si;
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity && modify_g){
+ /* find avg constitutive matrix from all vpts (change this later) */
+ get_rtf_at_vpts(E, m, lev, el, rtf2);
+ for(i=0;i<6;i++)
+ for(j=0;j<6;j++)
+ Duse[i][j]=0.0;
+ weight = 1./(2.*vpts);
+ for(i=1;i <= vpts;i++){ /* get vag const matrix */
+ off = (el-1)*vpts+i;
+ get_constitutive(Dtmp,lev,m,off,rtf2[1][i],rtf2[2][i],TRUE,E);
+ for(j=0;j<6;j++)
+ for(k=0;k<6;k++)
+ Duse[j][k] += Dtmp[j][k]*weight;
+ }
+ get_ba_p(&(E->N),&(E->GNX[lev][m][el]),&E->element_Cc, &E->element_Ccx,rtf,E->mesh.nsd,ba);
- for(a=1;a<=ends;a++) {
- for (i=1;i<=dims;i++)
- x[i]=E->GNX[lev][m][el].ppt[GNPXINDEX(2,a,1)]*E->element_Cc.ppt[BPINDEX(3,i,a,1)]
- + 2.0*ra*E->N.ppt[GNPINDEX(a,1)]*E->element_Cc.ppt[BPINDEX(3,i,a,1)]
- + ra*(E->GNX[lev][m][el].ppt[GNPXINDEX(0,a,1)]*E->element_Cc.ppt[BPINDEX(1,i,a,1)]
- +E->N.ppt[GNPINDEX(a,1)]*E->element_Ccx.ppt[BPXINDEX(1,i,1,a,1)]
- +ct*E->N.ppt[GNPINDEX(a,1)]*E->element_Cc.ppt[BPINDEX(1,i,a,1)]
- +si*(E->GNX[lev][m][el].ppt[GNPXINDEX(1,a,1)]*E->element_Cc.ppt[BPINDEX(2,i,a,1)]
- +E->N.ppt[GNPINDEX(a,1)]*E->element_Ccx.ppt[BPXINDEX(2,i,2,a,1)]));
+ /* assume single pressure point */
+ for(a = 1; a <= ends; a++){
+ for(i = 1; i <= dims; i++){
+ x[i] = 0.0;
+ for(k=0;k < 6;k++){
+ x[i] += Duse[0][k] * ba[a][1][i][k+1];
+ x[i] += Duse[1][k] * ba[a][1][i][k+1];
+ x[i] += Duse[2][k] * ba[a][1][i][k+1];
+ }
+ }
+ p=dims*(a-1);
+ elt_del[p ][0] = -x[1] * temp;
+ elt_del[p+1][0] = -x[2] * temp;
+ elt_del[p+2][0] = -x[3] * temp;
+
+ }
+ }else{
+#endif
+ ra = rtf[3][1];
+ si = 1.0/sin(rtf[1][1]);
+ ct = cos(rtf[1][1])*si;
- p=dims*(a-1);
- elt_del[p ][0] = -x[1] * temp;
- elt_del[p+1][0] = -x[2] * temp;
- elt_del[p+2][0] = -x[3] * temp;
-
+ /* old, isotropic part */
+ for(a=1;a<=ends;a++) {
+ for (i=1;i<=dims;i++)
+ x[i] = E->GNX[lev][m][el].ppt[GNPXINDEX(2,a,1)] * E->element_Cc.ppt[BPINDEX(3,i,a,1)] +
+ 2.0 * ra * E->N.ppt[GNPINDEX(a,1)]*E->element_Cc.ppt[BPINDEX(3,i,a,1)] +
+ ra *
+ (E->GNX[lev][m][el].ppt[GNPXINDEX(0,a,1)]*E->element_Cc.ppt[BPINDEX(1,i,a,1)] +
+ E->N.ppt[GNPINDEX(a,1)]*E->element_Ccx.ppt[BPXINDEX(1,i,1,a,1)] +
+ ct * E->N.ppt[GNPINDEX(a,1)] * E->element_Cc.ppt[BPINDEX(1,i,a,1)] +
+ si * (E->GNX[lev][m][el].ppt[GNPXINDEX(1,a,1)] * E->element_Cc.ppt[BPINDEX(2,i,a,1)] +
+ E->N.ppt[GNPINDEX(a,1)] * E->element_Ccx.ppt[BPXINDEX(2,i,2,a,1)]));
+ p=dims*(a-1);
+ elt_del[p ][0] = -x[1] * temp;
+ elt_del[p+1][0] = -x[2] * temp;
+ elt_del[p+2][0] = -x[3] * temp;
+
/* fprintf (E->fp,"B= %d %d %g %g %g %g %g\n",el,a,E->GDA[lev][m][el].ppt[1],E->GNX[lev][m][el].ppt[GNPXINDEX(0,a,1)],E->GNX[lev][m][el].ppt[GNPXINDEX(1,a,1)],elt_del[p][0],elt_del[p+1][0]);
- */
+ */
}
-
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ }
+#endif
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2010-09-22 23:48:53 UTC (rev 17211)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2010-09-24 03:35:36 UTC (rev 17212)
@@ -428,6 +428,8 @@
E->viscosity.zbase_layer[0] = E->viscosity.zbase_layer[1] = -999;
input_float_vector("z_layer",E->viscosity.num_mat,(E->viscosity.zbase_layer),m);
+
+
/* the start age and initial subduction history */
input_float("start_age",&(E->control.start_age),"0.0",m);
input_int("reset_startage",&(E->control.reset_startage),"0",m);
@@ -791,7 +793,10 @@
E->viscosity.zbase_layer[0] = E->viscosity.zlith;
E->viscosity.zbase_layer[1] = E->viscosity.z410;
E->viscosity.zbase_layer[2] = E->viscosity.zlm;
- E->viscosity.zbase_layer[3] = E->viscosity.zcmb;
+ E->viscosity.zbase_layer[3] = E->viscosity.zcmb; /* the lowest layers is never checked, really
+ if x3 < zlm, then the last layers gets assigned
+ i left this in for backward compatibility
+ */
}
return;
More information about the CIG-COMMITS
mailing list