[cig-commits] r7648 - in long/3D/Gale/trunk: . src/Gale/Utils/src
walter at geodynamics.org
walter at geodynamics.org
Thu Jul 12 11:11:51 PDT 2007
Author: walter
Date: 2007-07-12 11:11:51 -0700 (Thu, 12 Jul 2007)
New Revision: 7648
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.c
Log:
r1866 at earth: boo | 2007-07-12 11:06:33 -0700
Make friction work with distorted boundaries
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1861
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1866
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.c 2007-07-12 00:52:34 UTC (rev 7647)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.c 2007-07-12 18:11:51 UTC (rev 7648)
@@ -590,12 +590,12 @@
IndexSet *set = NULL;
Stream* warningStr = Journal_Register( Error_Type, self->type );
unsigned nDims;
- Grid* vertGrid;
+ Grid* grid;
nDims = Mesh_GetDimSize( self->_mesh );
- vertGrid = *(Grid**)ExtensionManager_Get( self->_mesh->info, self->_mesh,
- ExtensionManager_GetHandle( self->_mesh->info,
- "vertexGrid" ) );
+ grid = *(Grid**)ExtensionManager_Get( self->_mesh->info, self->_mesh,
+ ExtensionManager_GetHandle( self->_mesh->info,
+ "vertexGrid" ) );
/* We get the stress and pressure here. It requires that you
name your pressure and stress as PressureField and
StressField. I could not figure out how to get it to work
@@ -607,104 +607,316 @@
if(!self->pressure)
self->pressure = (FeVariable*)FieldVariable_Register_GetByName
( self->context->fieldVariable_Register, "PressureField" );
-
- switch (self->_wall) {
- case FrictionVC_Wall_Front:
- if ( nDims < 3 || !vertGrid->sizes[2] ) {
- set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
- }
- else {
- set = RegularMeshUtils_CreateGlobalFrontSet( self->_mesh );
- }
- break;
-
- case FrictionVC_Wall_Back:
- if ( nDims < 3 || !vertGrid->sizes[2] ) {
- set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
- }
- else {
- set = RegularMeshUtils_CreateGlobalBackSet( self->_mesh );
- }
- break;
-
- case FrictionVC_Wall_Top:
- if ( nDims < 2 || !vertGrid->sizes[1] ) {
- set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
- }
- else {
- set = RegularMeshUtils_CreateGlobalTopSet(self->_mesh);
- }
- break;
-
- case FrictionVC_Wall_Bottom:
- if ( nDims < 2 || !vertGrid->sizes[1] ) {
- set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
- }
- else {
- if(self->deviatoric_stress!=NULL && self->deviatoric_stress->dofLayout!=NULL)
+
+ /* If this is the first time this routine is called, the
+ stress may not have been initialized yet. In that case,
+ make everything stick. */
+ if(self->deviatoric_stress!=NULL && self->deviatoric_stress->dofLayout!=NULL)
+ {
+ Node_LocalIndex n_i;
+ unsigned nNodes, direction, boundary, basis[2];
+ IJK ijk;
+
+ /* Note that we select the order of the basis vectors to point inwards */
+ switch (self->_wall) {
+ case FrictionVC_Wall_Front:
+ direction=K_AXIS;
+ boundary=grid->sizes[direction]-1;
+ basis[0]=J_AXIS;
+ basis[1]=I_AXIS;
+ break;
+ case FrictionVC_Wall_Back:
+ direction=K_AXIS;
+ boundary=0;
+ basis[0]=I_AXIS;
+ basis[1]=J_AXIS;
+ break;
+ case FrictionVC_Wall_Top:
+ direction=J_AXIS;
+ boundary=grid->sizes[direction]-1;
+ basis[0]=I_AXIS;
+ basis[1]=K_AXIS;
+ break;
+ case FrictionVC_Wall_Bottom:
+ direction=J_AXIS;
+ boundary=0;
+ basis[0]=K_AXIS;
+ basis[1]=I_AXIS;
+ break;
+ case FrictionVC_Wall_Left:
+ direction=I_AXIS;
+ boundary=0;
+ basis[0]=K_AXIS;
+ basis[1]=J_AXIS;
+ break;
+ case FrictionVC_Wall_Right:
+ direction=I_AXIS;
+ boundary=grid->sizes[direction]-1;
+ basis[0]=J_AXIS;
+ basis[1]=K_AXIS;
+ break;
+ case FrictionVC_Wall_Size:
+ default:
+ assert(0);
+ break;
+ }
+
+ nNodes = Mesh_GetDomainSize( self->_mesh, MT_VERTEX );
+ set = IndexSet_New( nNodes );
+ for( n_i = 0; n_i < nNodes; n_i++ ) {
+ RegularMeshUtils_Node_1DTo3D
+ ( self->_mesh, Mesh_DomainToGlobal
+ ( self->_mesh, MT_VERTEX, n_i ), ijk );
+ if(ijk[direction]==boundary)
+ {
+ double str[6], p, normal_stress, tangential_norm;
+ FeVariable_GetValueAtNode(self->deviatoric_stress,n_i,str);
+ FeVariable_GetValueAtNode(self->pressure,n_i,&p);
+
+ if(nDims==2)
{
- Node_LocalIndex n_i;
- unsigned nNodes;
- IJK ijk;
+ double x1[2],x2[2],surface[2],n[2],tangential_stress[2];
+ int sign;
+ unsigned n_temp, ijk_temp[3];
+ 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)
+ {
+ basis[0]=basis[1];
+ sign=-1;
+ }
+ /* Get x1 */
+ if(ijk[basis[0]]==0)
+ {
+ x1[0]=self->_mesh->verts[n_i][0];
+ x1[1]=self->_mesh->verts[n_i][1];
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[basis[0]]-=1;
+ n_temp=
+ RegularMeshUtils_Node_3DTo1D(self->_mesh,ijk_temp);
+ x1[0]=self->_mesh->verts[n_temp][0];
+ x1[1]=self->_mesh->verts[n_temp][1];
+ }
+ /* Get x2 */
+ if(ijk[basis[0]]==grid->sizes[basis[0]]-1)
+ {
+ x2[0]=self->_mesh->verts[n_i][0];
+ x2[1]=self->_mesh->verts[n_i][1];
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[basis[0]]+=1;
+ n_temp=
+ RegularMeshUtils_Node_3DTo1D(self->_mesh,ijk_temp);
+ x2[0]=self->_mesh->verts[n_temp][0];
+ x2[1]=self->_mesh->verts[n_temp][1];
+ }
+ /* Get the surface */
+ Vec_Sub2D(surface,x2,x1);
+ n[0]=sign*surface[1];
+ n[1]=sign*surface[0];
+ Vec_Norm2D(n,n);
- nNodes = Mesh_GetDomainSize( self->_mesh, MT_VERTEX );
- set = IndexSet_New( nNodes );
- for( n_i = 0; n_i < nNodes; n_i++ ) {
- RegularMeshUtils_Node_1DTo3D( self->_mesh, Mesh_DomainToGlobal( self->_mesh, MT_VERTEX, n_i ), ijk );
- if(ijk[1]==0)
- {
- double str[6], p;
- FeVariable_GetValueAtNode(self->deviatoric_stress,n_i,str);
- FeVariable_GetValueAtNode(self->pressure,n_i,&p);
+ /* Normal stress = n . stress . n */
- /* For now, assume that the bottom is flat
- and not inclined.*/
+ /* The ordering for components of stress is: xx,
+ yy, xy */
- /* The ordering for components of stress
- is: xx, yy, zz, xy, xz, zz */
+ normal_stress=n[0]*n[0]*str[0]
+ + n[1]*n[1]*str[1]
+ + 2*n[1]*n[0]*str[2];
- if((nDims==2
- && (str[2]==0 ||
- self->friction*(p-str[1])>fabs(str[2])))
- || (nDims==3
- && ((str[3]==0 && str[4]==0)
- || self->friction*(p-str[1])>sqrt(str[3]*str[3] + str[4]*str[4]))))
- IndexSet_Add(set,n_i);
- }
- }
+ /* Tangential stress = sigma . n - n * normal_stress
+ Note that it is a vector */
+
+ tangential_stress[0]=str[0]*n[0]
+ + str[2]*n[1]
+ - n[0]*normal_stress;
+ tangential_stress[1]=str[1]*n[1]
+ + str[2]*n[0]
+ - n[1]*normal_stress;
+
+ tangential_norm=Vec_Mag2D(tangential_stress);
}
else
{
- set = RegularMeshUtils_CreateGlobalBottomSet(self->_mesh);
+ double x1[3],x2[3],x3[3],x4[3],line1[3],line2[3],n[3],
+ tangential_stress[3];
+ unsigned n_temp, ijk_temp[3];
+
+ /* Get x1 */
+ if(ijk[basis[0]]==0)
+ {
+ x1[0]=self->_mesh->verts[n_i][0];
+ x1[1]=self->_mesh->verts[n_i][1];
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[basis[0]]-=1;
+ n_temp=
+ RegularMeshUtils_Node_3DTo1D(self->_mesh,ijk_temp);
+ x1[0]=self->_mesh->verts[n_temp][0];
+ x1[1]=self->_mesh->verts[n_temp][1];
+ }
+ /* Get x2 */
+ if(ijk[basis[0]]==grid->sizes[basis[0]]-1)
+ {
+ x2[0]=self->_mesh->verts[n_i][0];
+ x2[1]=self->_mesh->verts[n_i][1];
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[basis[0]]+=1;
+ n_temp=
+ RegularMeshUtils_Node_3DTo1D(self->_mesh,ijk_temp);
+ x2[0]=self->_mesh->verts[n_temp][0];
+ x2[1]=self->_mesh->verts[n_temp][1];
+ }
+
+ /* Get x3 */
+ if(ijk[basis[1]]==0)
+ {
+ x3[0]=self->_mesh->verts[n_i][0];
+ x3[1]=self->_mesh->verts[n_i][1];
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[basis[1]]-=1;
+ n_temp=
+ RegularMeshUtils_Node_3DTo1D(self->_mesh,ijk_temp);
+ x3[0]=self->_mesh->verts[n_temp][0];
+ x3[1]=self->_mesh->verts[n_temp][1];
+ }
+ /* Get x4 */
+ if(ijk[basis[1]]==grid->sizes[basis[1]]-1)
+ {
+ x4[0]=self->_mesh->verts[n_i][0];
+ x4[1]=self->_mesh->verts[n_i][1];
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[basis[1]]+=1;
+ n_temp=
+ RegularMeshUtils_Node_3DTo1D(self->_mesh,ijk_temp);
+ x4[0]=self->_mesh->verts[n_temp][0];
+ x4[1]=self->_mesh->verts[n_temp][1];
+ }
+
+ /* Compute the normal. */
+ Vec_Sub2D(line1,x2,x1);
+ Vec_Sub2D(line2,x4,x3);
+ Vec_Cross3D(n,line1,line2);
+ Vec_Norm3D(n,n);
+
+ /* Normal stress = n . stress . n */
+
+ /* The ordering for components of stress is: xx,
+ yy, zz, xy, xz, zz */
+
+ normal_stress=n[0]*n[0]*str[0]
+ + n[1]*n[1]*str[1]
+ + n[2]*n[2]*str[2]
+ + 2*n[1]*n[0]*str[3]
+ + 2*n[2]*n[0]*str[4]
+ + 2*n[2]*n[1]*str[5];
+
+ /* Tangential stress = sigma . n - n * normal_stress
+ Note that it is a vector */
+
+ tangential_stress[0]=str[0]*n[0]
+ + str[3]*n[1] + str[4]*n[2]
+ - n[0]*normal_stress;
+ tangential_stress[1]=str[3]*n[0]
+ + str[1]*n[1] + str[5]*n[2]
+ - n[1]*normal_stress;
+ tangential_stress[2]=str[4]*n[0]
+ + str[5]*n[1] + str[2]*n[2]
+ - n[2]*normal_stress;
+
+ tangential_norm=Vec_Mag3D(tangential_stress);
}
- }
- break;
-
- case FrictionVC_Wall_Left:
- if ( nDims < 1 ) {
- set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
- }
- else {
- set = RegularMeshUtils_CreateGlobalLeftSet(self->_mesh);
- }
- break;
-
- case FrictionVC_Wall_Right:
- if( nDims < 1 ) {
- set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
- }
- else {
- set = RegularMeshUtils_CreateGlobalRightSet(self->_mesh);
- }
- break;
-
- case FrictionVC_Wall_Size:
- default:
- assert(0);
- break;
- }
-
- return set;
+
+ /* Finally, add in the point if the frictional force
+ keeps the material pinned to the boundary. */
+ if(self->friction*(p-normal_stress)>=tangential_norm)
+ IndexSet_Add(set,n_i);
+ }
+ }
+ }
+ else
+ {
+ switch (self->_wall) {
+ case FrictionVC_Wall_Front:
+ if ( nDims < 3 || !grid->sizes[2] ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
+ }
+ else {
+ set = RegularMeshUtils_CreateGlobalFrontSet( self->_mesh );
+ }
+ break;
+
+ case FrictionVC_Wall_Back:
+ if ( nDims < 3 || !grid->sizes[2] ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
+ }
+ else {
+ set = RegularMeshUtils_CreateGlobalBackSet( self->_mesh );
+ }
+ break;
+
+ case FrictionVC_Wall_Top:
+ if ( nDims < 2 || !grid->sizes[1] ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
+ }
+ else {
+ set = RegularMeshUtils_CreateGlobalTopSet(self->_mesh);
+ }
+ break;
+
+ case FrictionVC_Wall_Bottom:
+ if ( nDims < 2 || !grid->sizes[1] ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
+ }
+ else {
+ set = RegularMeshUtils_CreateGlobalBottomSet(self->_mesh);
+ }
+ break;
+
+ case FrictionVC_Wall_Left:
+ if ( nDims < 1 ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
+ }
+ else {
+ set = RegularMeshUtils_CreateGlobalLeftSet(self->_mesh);
+ }
+ break;
+
+ case FrictionVC_Wall_Right:
+ if( nDims < 1 ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
+ }
+ else {
+ set = RegularMeshUtils_CreateGlobalRightSet(self->_mesh);
+ }
+ break;
+
+ case FrictionVC_Wall_Size:
+ default:
+ assert(0);
+ break;
+ }
+ }
+ return set;
}
More information about the cig-commits
mailing list