[cig-commits] r11974 - in long/3D/Gale/trunk: . src/Gale/Utils/src
walter at geodynamics.org
walter at geodynamics.org
Thu May 15 12:07:05 PDT 2008
Author: walter
Date: 2008-05-15 12:07:05 -0700 (Thu, 15 May 2008)
New Revision: 11974
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c
long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c
Log:
r2154 at earth: boo | 2008-05-15 09:32:59 -0700
Go back to having the kinetic friction inside the nonlinear loop
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2153
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2154
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c 2008-05-15 19:07:00 UTC (rev 11973)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c 2008-05-15 19:07:05 UTC (rev 11974)
@@ -377,16 +377,196 @@
if(ijk[direction]==boundary)
{
+ unsigned overcount=StressBC_get_overcount(dim,ijk,grid->sizes);
+ double v[3], v_surface[3], v_surface_diff[3], normal_stress,
+ tangential_stress[3], tangential_norm, friction_coefficient;
+ Bool include_boundary[3];
+
+ double p;
+
double *force_vector=StaticFrictionVC_get_force_vector();
+ Stream* errorStr = Journal_Register( Error_Type, self->type );
+ switch(self->friction_entry.type)
+ {
+ case StressBC_Double:
+ friction_coefficient=self->friction_entry.DoubleValue;
+ break;
+ case StressBC_ConditionFunction:
+ cf=self->conFunc_Register->
+ _cf[self->friction_entry.CFIndex];
+
+ /* We use a variable number of zero "0", because we don't
+ use the variable number and that one is always going to
+ exist. */
+ ConditionFunction_Apply(cf,elementNodes[eNode_I],
+ 0,self->context,&friction_coefficient);
+ break;
+ }
+
for(d=0;d<dim;++d)
{
- elForceVec[ eNode_I*nodeDofCount + d]+=
- force_vector[elementNodes[eNode_I]*dim+d];
+ switch(self->v_entry[d].type)
+ {
+ case StressBC_Double:
+ v_boundary[d]=self->v_entry[d].DoubleValue;
+ break;
+ case StressBC_ConditionFunction:
+ cf=self->conFunc_Register
+ ->_cf[self->v_entry[d].CFIndex];
+ ConditionFunction_Apply(cf,elementNodes[eNode_I],
+ 0,self->context,v_boundary+d);
+ break;
+ }
+ }
+ deviatoric_stress =
+ (FeVariable*)FieldVariable_Register_GetByName
+ ( self->context->fieldVariable_Register,
+ self->deviatoric_stress_name );
+ pressure = (FeVariable*)FieldVariable_Register_GetByName
+ ( self->context->fieldVariable_Register,
+ self->pressure_name );
+
+ Journal_Firewall( ( pressure!=NULL ), errorStr,
+ "Error: In StressVC, the name provided for the pressure field \"%s\" does not exist\n",
+ self->pressure_name);
+ Journal_Firewall( ( deviatoric_stress!=NULL ), errorStr,
+ "Error: In StressVC, the name provided for the stress field \"%s\" does not exist\n",
+ self->deviatoric_stress_name);
+
+ /* We always include the boundaries when calculating norms and
+ stresses. The include_boundary argument is really for
+ StaticFrictionVC's. */
+ for(d=0;d<dim;d++)
+ include_boundary[d]=True;
+ StaticFrictionVC_get_interface_stresses
+ (deviatoric_stress,mesh,grid->sizes,elementNodes[eNode_I],
+ ijk,dim,basis,include_boundary,
+ include_boundary,&normal_stress,tangential_stress,&tangential_norm,n);
- printf("kinetic %d %d %d %g\n",elementNodes[eNode_I],eNode_I,d,
- force_vector[elementNodes[eNode_I]*dim+d]);
+ {
+/* double p; */
+ FeVariable_GetValueAtNode(pressure,lElement_I,&p);
+ normal_stress+=p;
+ }
+
+ velocity = (FeVariable*)FieldVariable_Register_GetByName
+ ( self->context->fieldVariable_Register,
+ self->velocity_name );
+ Journal_Firewall( ( velocity!=NULL ), errorStr,
+ "Error: In KineticFriction, the name provided for the velocity field \"%s\" does not exist\n",
+ self->velocity_name);
+ FeVariable_GetValueAtNode(velocity,elementNodes[eNode_I],v);
+
+ /* Get the magnitude of the velocity parallel to the surface. */
+ if(dim==2)
+ {
+ Vec_Sub2D(v_diff,v_boundary,v);
+ v_dot_n=Vec_Dot2D(v_diff,n);
+ Vec_Scale2D(v_perpendicular,n,v_dot_n);
+ Vec_Sub2D(v_surface_diff,v_diff,v_perpendicular);
+ v_norm=Vec_Mag2D(v_surface_diff);
}
+ else
+ {
+ Vec_Sub3D(v_diff,v_boundary,v);
+ v_dot_n=Vec_Dot3D(v,n);
+ Vec_Scale3D(v_perpendicular,n,v_dot_n);
+ Vec_Sub3D(v_surface_diff,v,v_perpendicular);
+ v_norm=Vec_Mag3D(v_surface_diff);
+ }
+
+ StaticFrictionVC_get_force_vector(force_vector);
+
+/* if(force_vector[elementNodes[eNode_I]]>0 */
+/* && v_norm>1.0e-11) */
+/* { */
+/* for(d=0;d<dim;++d) */
+/* { */
+/* elForceVec[ eNode_I*nodeDofCount + d]+=force_vector[elementNodes[eNode_I]]*(v_surface_diff[d]/v_norm)*area/overcount; */
+/* } */
+/* printf("kinetic slip %d %d %d %g %g %g %g\n", */
+/* lElement_I,eNode_I,elementNodes[eNode_I], */
+/* force_vector[elementNodes[eNode_I]]*area/overcount, */
+/* v_norm,v_surface_diff[0],v_surface_diff[1]); */
+/* } */
+/* else */
+/* { */
+/* printf("kinetic stick %d %d %d %g %g %g %g\n", */
+/* lElement_I,eNode_I,elementNodes[eNode_I], */
+/* force_vector[elementNodes[eNode_I]]*area/overcount, */
+/* v_norm,v_surface_diff[0],v_surface_diff[1]); */
+/* } */
+
+ friction_force=friction_coefficient*normal_stress*area/overcount;
+
+/* if(!(friction_force<0 || v_norm==0)) */
+ if(friction_force>0
+/* && friction_coefficient*normal_stress<tangential_norm */
+ && v_norm>1.0e-11)
+ {
+ for(d=0;d<dim;++d)
+ {
+ double alpha=0.5;
+ double old_force=
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d];
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d]=
+ friction_force*v_surface_diff[d]/v_norm;
+ if(old_force==0)
+ alpha=1.0;
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d]=
+ (friction_force*v_surface_diff[d]/v_norm)*alpha
+ + old_force*(1-alpha);
+
+ if(force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d]*old_force<0)
+ {
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d]*=1.0e-20;
+ printf("kinetic reduced\n");
+ }
+/* force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d]= */
+/* (friction_force*v_surface_diff[d]/v_norm-old_force) */
+/* - tangential_stress[d]*area/overcount; */
+
+ elForceVec[ eNode_I*nodeDofCount + d]+=
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d];
+
+/* elForceVec[ eNode_I*nodeDofCount + d]+=friction_force*v_surface_diff[d]/v_norm; */
+/* elForceVec[ eNode_I*nodeDofCount + d]+=friction_force*v_surface_diff[d]/v_norm - tangential_stress[d]*area/overcount; */
+/* node_force[ eNode_I*nodeDofCount + d]+=friction_force*v_surface_diff[d]/v_norm; */
+/* node_force[ eNode_I*nodeDofCount + d]+=friction_force*v_surface_diff[d]/v_norm - tangential_stress[d]*area/overcount; */
+/* for(d=0;d<dim;++d) */
+
+ printf("kinetic slip %d %d %d %d %g %g %g %g %g %g %g %g %d %g %g\n",
+ (lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d,
+ lElement_I,eNode_I,d,friction_force,old_force,
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d],
+ v_norm,v_surface_diff[0],v_surface_diff[1],
+ normal_stress,p,overcount,tangential_stress[0],
+ tangential_stress[1]);
+ }
+ }
+ else
+ {
+ for(d=0;d<dim;++d)
+ {
+ double old_force=
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d];
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d]=0;
+
+/* force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d]= */
+/* - (tangential_stress[d]*area/overcount-old_force); */
+/* elForceVec[ eNode_I*nodeDofCount + d]+= */
+/* force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d]; */
+
+ printf("kinetic stick %d %d %d %d %g %g %g %g %g %g %g %g %d %g %g\n",
+ (lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d,
+ lElement_I,eNode_I,d,friction_force,old_force,
+ force_vector[(lElement_I*elementNodeCount+eNode_I)*nodeDofCount+d],
+ v_norm,v_surface_diff[0],v_surface_diff[1],
+ normal_stress,p,overcount,tangential_stress[0],
+ tangential_stress[1]);
+ }
+ }
}
}
}
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c 2008-05-15 19:07:00 UTC (rev 11973)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c 2008-05-15 19:07:05 UTC (rev 11974)
@@ -617,7 +617,7 @@
Stream* errorStr = Journal_Register( Error_Type, self->type );
unsigned nDims;
Grid* grid;
- FeVariable *deviatoric_stress, *pressure, *velocity;
+ FeVariable *deviatoric_stress, *pressure;
unsigned nnn;
nDims = Mesh_GetDimSize( self->_mesh );
@@ -629,18 +629,6 @@
pressure = (FeVariable*)FieldVariable_Register_GetByName
( self->context->fieldVariable_Register, self->pressure_name );
- velocity = (FeVariable*)FieldVariable_Register_GetByName
- ( self->context->fieldVariable_Register,
- "VelocityField" );
-
-/* velocity = (FeVariable*)FieldVariable_Register_GetByName */
-/* ( self->context->fieldVariable_Register, */
-/* self->velocity_name ); */
-
- Journal_Firewall( ( velocity!=NULL ), errorStr,
- "Error: In StaticFriction, the name provided for the velocity field \"%s\" does not exist\n",
- "VelocityField");
-
Journal_Firewall( ( pressure!=NULL ), errorStr,
"Error: In StaticFrictionVC, the name provided for the pressure field \"%s\" does not exist\n",
self->pressure_name);
@@ -715,8 +703,6 @@
}
else
printf("static slip initial %d\n",n_i);
-
-/* force_vector[n_i*nDims+0]=-0.5*0.5/32; */
}
/* Finally, add in the point if the frictional force
keeps the material pinned to the boundary. */
@@ -726,47 +712,16 @@
self->include_lower_boundary,
self->include_upper_boundary,
&normal_stress,tangential_stress,&tangential_norm,n)
-/* && 0.50*(p+normal_stress)>=tangential_norm) */
- && 0.95*self->friction*(p+normal_stress)>=tangential_norm)
+ && self->friction*(p+normal_stress)>=tangential_norm)
{
- int d;
IndexSet_Add(set,n_i);
- printf("static stick %d %g %g %g\n",
+ printf("static stick %d %lf %lf %lf\n",
n_i,p,normal_stress,tangential_norm);
- for(d=0;d<nDims;++d)
- force_vector[n_i*nDims+d]=0.0;
}
else
{
- double v_boundary;
- double v[3];
- double *coord;
- double area=4/64.0;
- int overcount=2;
-
- coord = Mesh_GetVertex(self->_mesh,n_i);
-
-/* v_boundary=0; */
-
- if(coord[0]<0.1)
- v_boundary=0;
- else
- v_boundary=6.94444444444e-6;
-
- FeVariable_GetValueAtNode(velocity,n_i,v);
-
- if(fabs(v_boundary-v[0])>1.0e-11 && p+normal_stress>0)
- {
- if(v_boundary>v[0])
- force_vector[n_i*nDims+0]=self->friction*(p+normal_stress)*area/overcount;
- else
- force_vector[n_i*nDims+0]=-self->friction*(p+normal_stress)*area/overcount;
- }
- else
- force_vector[n_i*nDims+0]=0;
-
- printf("static slip %d %g %g %g %g %g %g %g %g %g\n",
- n_i,coord[0],p,normal_stress,tangential_norm,v_boundary,v[0],v[1],force_vector[n_i*nDims+0],p+normal_stress);
+ printf("static slip %d %lf %lf %lf\n",
+ n_i,p,normal_stress,tangential_norm);
}
}
}
@@ -1067,7 +1022,7 @@
*tangential_norm=Vec_Mag2D(tangential_stress);
-/* printf("interface %g %g %g %g %d\n", */
+/* printf("interface %lf %lf %lf %lf %d\n", */
/* n[0],n[1],surface[0],surface[1],sign); */
}
More information about the cig-commits
mailing list