[cig-commits] r11973 - in long/3D/Gale/trunk: . src/Gale/Utils/src
walter at geodynamics.org
walter at geodynamics.org
Thu May 15 12:07:00 PDT 2008
Author: walter
Date: 2008-05-15 12:07:00 -0700 (Thu, 15 May 2008)
New Revision: 11973
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:
r2153 at earth: boo | 2008-05-15 09:30:54 -0700
Awful hack to take kinetic friction out of the nonlinear loop. Does not work with the extension benchmark
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2146
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2153
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c 2008-05-15 06:00:17 UTC (rev 11972)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c 2008-05-15 19:07:00 UTC (rev 11973)
@@ -377,196 +377,16 @@
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)
{
- 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);
+ elForceVec[ eNode_I*nodeDofCount + 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);
+ printf("kinetic %d %d %d %g\n",elementNodes[eNode_I],eNode_I,d,
+ force_vector[elementNodes[eNode_I]*dim+d]);
}
- 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 06:00:17 UTC (rev 11972)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c 2008-05-15 19:07:00 UTC (rev 11973)
@@ -617,7 +617,7 @@
Stream* errorStr = Journal_Register( Error_Type, self->type );
unsigned nDims;
Grid* grid;
- FeVariable *deviatoric_stress, *pressure;
+ FeVariable *deviatoric_stress, *pressure, *velocity;
unsigned nnn;
nDims = Mesh_GetDimSize( self->_mesh );
@@ -629,6 +629,18 @@
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);
@@ -703,6 +715,8 @@
}
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. */
@@ -712,16 +726,47 @@
self->include_lower_boundary,
self->include_upper_boundary,
&normal_stress,tangential_stress,&tangential_norm,n)
- && self->friction*(p+normal_stress)>=tangential_norm)
+/* && 0.50*(p+normal_stress)>=tangential_norm) */
+ && 0.95*self->friction*(p+normal_stress)>=tangential_norm)
{
+ int d;
IndexSet_Add(set,n_i);
- printf("static stick %d %lf %lf %lf\n",
+ printf("static stick %d %g %g %g\n",
n_i,p,normal_stress,tangential_norm);
+ for(d=0;d<nDims;++d)
+ force_vector[n_i*nDims+d]=0.0;
}
else
{
- printf("static slip %d %lf %lf %lf\n",
- n_i,p,normal_stress,tangential_norm);
+ 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);
}
}
}
@@ -1022,7 +1067,7 @@
*tangential_norm=Vec_Mag2D(tangential_stress);
-/* printf("interface %lf %lf %lf %lf %d\n", */
+/* printf("interface %g %g %g %g %d\n", */
/* n[0],n[1],surface[0],surface[1],sign); */
}
More information about the cig-commits
mailing list