[cig-commits] r11967 - in long/3D/Gale/trunk: . src/Gale/Utils/src
walter at geodynamics.org
walter at geodynamics.org
Wed May 14 22:59:46 PDT 2008
Author: walter
Date: 2008-05-14 22:59:46 -0700 (Wed, 14 May 2008)
New Revision: 11967
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
long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.h
Log:
r2141 at earth: boo | 2008-05-13 01:28:27 -0700
Add the ability to specify whether static friction sticks or slips on the first time step. Make StaticFrictionVC_get_interface_stresses return tangential stress, including fixing a bug with the direction of the stress. Some hacks to KineticFriction that do not seem to work.
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2138
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2141
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c 2008-05-14 16:23:22 UTC (rev 11966)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c 2008-05-15 05:59:46 UTC (rev 11967)
@@ -342,7 +342,15 @@
unsigned direction, boundary, basis[2], d;
FeVariable *pressure, *deviatoric_stress, *velocity;
Node_DomainIndex* elementNodes = NULL;
-
+
+ /* Hard code the dimension and number of nodes */
+ double node_force[3*4];
+ Bool apply_force;
+
+ apply_force=True;
+ for(eNode_I=0;eNode_I<12;++eNode_I)
+ node_force[eNode_I]=0;
+
pressure=deviatoric_stress=velocity=NULL;
elementType = FeMesh_GetElementType( mesh, lElement_I );
@@ -371,8 +379,11 @@
{
unsigned overcount=StressBC_get_overcount(dim,ijk,grid->sizes);
double v[3], v_surface[3], v_surface_diff[3], normal_stress,
- tangential_norm, friction_coefficient;
+ tangential_stress[3], tangential_norm, friction_coefficient;
Bool include_boundary[3];
+
+ double p;
+
Stream* errorStr = Journal_Register( Error_Type, self->type );
switch(self->friction_entry.type)
{
@@ -429,10 +440,10 @@
StaticFrictionVC_get_interface_stresses
(deviatoric_stress,mesh,grid->sizes,elementNodes[eNode_I],
ijk,dim,basis,include_boundary,
- include_boundary,&normal_stress,&tangential_norm,n);
+ include_boundary,&normal_stress,tangential_stress,&tangential_norm,n);
{
- double p;
+/* double p; */
FeVariable_GetValueAtNode(pressure,lElement_I,&p);
normal_stress+=p;
}
@@ -464,10 +475,44 @@
}
friction_force=friction_coefficient*normal_stress*area/overcount;
- if(!(friction_force<0 || v_norm==0))
+/* 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)
+ {
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 %g %g %g %g %g %g %d %g %g\n",
+ lElement_I,eNode_I,friction_force,
+ v_norm,v_surface_diff[0],v_surface_diff[1],
+ normal_stress,p,overcount,tangential_stress[0],
+ tangential_stress[1]);
+ }
+ else
+ {
+ printf("kinetic stick %d %d %g %g %g %g %g %g %d %g %g\n",
+ lElement_I,eNode_I,friction_force,
+ v_norm,v_surface_diff[0],v_surface_diff[1],
+ normal_stress,p,overcount,tangential_stress[0],
+ tangential_stress[1]);
+ apply_force=False;
+ }
}
}
+
+/* if(apply_force) */
+/* { */
+/* for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) { */
+/* for(d=0;d<dim;++d) */
+/* elForceVec[ eNode_I*nodeDofCount + d]+= */
+/* node_force[eNode_I*nodeDofCount + d]; */
+/* } */
+/* } */
}
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c 2008-05-14 16:23:22 UTC (rev 11966)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c 2008-05-15 05:59:46 UTC (rev 11967)
@@ -262,6 +262,7 @@
self->include_lower_boundary[i]=True;
self->include_upper_boundary[i]=True;
}
+ self->stick_on_first_step=True;
self->friction=0;
}
@@ -323,6 +324,7 @@
self->include_upper_boundary[I_AXIS]=Dictionary_Entry_Value_AsBool_with_Default(Dictionary_Entry_Value_GetMember(vcDictVal, "includeUpperX"),True);
self->include_upper_boundary[J_AXIS]=Dictionary_Entry_Value_AsBool_with_Default(Dictionary_Entry_Value_GetMember(vcDictVal, "includeUpperY"),True);
self->include_upper_boundary[K_AXIS]=Dictionary_Entry_Value_AsBool_with_Default(Dictionary_Entry_Value_GetMember(vcDictVal, "includeUpperZ"),True);
+ self->stick_on_first_step=Dictionary_Entry_Value_AsBool_with_Default(Dictionary_Entry_Value_GetMember(vcDictVal, "StickOnFirstStep"),True);
self->friction = Dictionary_Entry_Value_AsDouble(Dictionary_Entry_Value_GetMember(vcDictVal, "StaticFriction"));
for (entry_I = 0; entry_I < self->_entryCount; entry_I++)
@@ -423,6 +425,7 @@
self->include_lower_boundary[i]=True;
self->include_upper_boundary[i]=True;
}
+ self->stick_on_first_step=True;
self->friction=0;
}
}
@@ -548,6 +551,7 @@
newStaticFrictionVC->include_lower_boundary[i]=self->include_lower_boundary[i];
newStaticFrictionVC->include_upper_boundary[i]=self->include_upper_boundary[i];
}
+ newStaticFrictionVC->stick_on_first_step = self->stick_on_first_step;
newStaticFrictionVC->friction = self->friction;
if( deep ) {
@@ -651,7 +655,8 @@
( self->_mesh, MT_VERTEX, n_i ), ijk );
if(ijk[direction]==boundary)
{
- double normal_stress, tangential_norm, n[3], p, temp;
+ double normal_stress, tangential_stress[3],
+ tangential_norm, n[3], p, temp;
/* Average the pressure over neighboring cells to
get the pressure at the node. */
@@ -667,19 +672,37 @@
}
p/=nInc;
+ /* If the pressure is zero, then we are at the first
+ time step. So we don't compute interface
+ stresses and just use the supplied initial
+ guess */
+ if(p==0)
+ {
+ if(self->stick_on_first_step)
+ {
+ IndexSet_Add(set,n_i);
+ printf("static stick initial %d\n",n_i);
+ }
+ else
+ printf("static slip initial %d\n",n_i);
+ }
/* Finally, add in the point if the frictional force
keeps the material pinned to the boundary. */
- if(StaticFrictionVC_get_interface_stresses
- (deviatoric_stress,self->_mesh,
- grid->sizes, n_i, ijk, nDims, basis,
- self->include_lower_boundary,
- self->include_upper_boundary,
- &normal_stress,&tangential_norm,n)
- && (self->friction*(p+normal_stress)>=tangential_norm
- || p+normal_stress==0))
+ else if(StaticFrictionVC_get_interface_stresses
+ (deviatoric_stress,self->_mesh,
+ grid->sizes, n_i, ijk, nDims, basis,
+ self->include_lower_boundary,
+ self->include_upper_boundary,
+ &normal_stress,tangential_stress,&tangential_norm,n)
+ && self->friction*(p+normal_stress)>=tangential_norm)
{
IndexSet_Add(set,n_i);
+ printf("static stick %d %lf %lf %lf\n",
+ n_i,p,normal_stress,tangential_norm);
}
+ else
+ printf("static slip %d %lf %lf %lf\n",
+ n_i,p,normal_stress,tangential_norm);
}
}
}
@@ -886,6 +909,7 @@
Bool include_lower_boundary[],
Bool include_upper_boundary[],
double *normal_stress,
+ double tangential_stress[],
double *tangential_norm,
double n[])
{
@@ -899,22 +923,24 @@
if(nDims==2)
{
- double x1[2],x2[2],surface[2],tangential_stress[2];
+ double x1[2],x2[2],surface[2];
int sign;
- unsigned n_temp, ijk_temp[3];
+ unsigned n_temp, ijk_temp[3], new_basis[2];
sign=1;
/* If the first direction is z, then use the
second direction and reverse it. This is
only for 2D. */
- if(basis[0]==K_AXIS)
+ new_basis[0]=basis[0];
+ new_basis[1]=basis[1];
+ if(new_basis[0]==K_AXIS)
{
- basis[0]=basis[1];
+ new_basis[0]=new_basis[1];
sign=-1;
}
/* Get x1 */
- if(ijk[basis[0]]==0)
+ if(ijk[new_basis[0]]==0)
{
- if(include_lower_boundary[basis[0]]==True)
+ if(include_lower_boundary[new_basis[0]]==True)
{
x1[0]=mesh->verts[n_i][0];
x1[1]=mesh->verts[n_i][1];
@@ -925,15 +951,15 @@
else
{
Vec_Set2D(ijk_temp,ijk);
- ijk_temp[basis[0]]-=1;
+ ijk_temp[new_basis[0]]-=1;
n_temp=RegularMeshUtils_Node_3DTo1D(mesh,ijk_temp);
x1[0]=mesh->verts[n_temp][0];
x1[1]=mesh->verts[n_temp][1];
}
/* Get x2 */
- if(ijk[basis[0]]==sizes[basis[0]]-1)
+ if(ijk[new_basis[0]]==sizes[new_basis[0]]-1)
{
- if(include_upper_boundary[basis[0]]==True)
+ if(include_upper_boundary[new_basis[0]]==True)
{
x2[0]=mesh->verts[n_i][0];
x2[1]=mesh->verts[n_i][1];
@@ -944,7 +970,7 @@
else
{
Vec_Set2D(ijk_temp,ijk);
- ijk_temp[basis[0]]+=1;
+ ijk_temp[new_basis[0]]+=1;
n_temp=RegularMeshUtils_Node_3DTo1D(mesh,ijk_temp);
x2[0]=mesh->verts[n_temp][0];
x2[1]=mesh->verts[n_temp][1];
@@ -976,10 +1002,13 @@
*tangential_norm=Vec_Mag2D(tangential_stress);
+/* printf("interface %lf %lf %lf %lf %d\n", */
+/* n[0],n[1],surface[0],surface[1],sign); */
+
}
else
{
- double x1[3],x2[3],x3[3],x4[3],line1[3],line2[3],tangential_stress[3];
+ double x1[3],x2[3],x3[3],x4[3],line1[3],line2[3];
unsigned n_temp, ijk_temp[3];
/* Get x1 */
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.h 2008-05-14 16:23:22 UTC (rev 11966)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.h 2008-05-15 05:59:46 UTC (rev 11967)
@@ -64,6 +64,7 @@
char *deviatoric_stress_name, *pressure_name; \
Bool include_lower_boundary[3]; \
Bool include_upper_boundary[3]; \
+ Bool stick_on_first_step; \
double friction;
struct _StaticFrictionVC { __StaticFrictionVC };
@@ -206,6 +207,7 @@
Bool include_lower_boundary[],
Bool include_upper_boundary[],
double *normal_stress,
+ double tangential_stress[],
double *tangential_norm,
double n[]);
More information about the cig-commits
mailing list