[cig-commits] r7696 - in long/3D/Gale/trunk: . src/Gale/Utils/src
src/PICellerator/Utils/src src/StgFEM/SLE/SystemSetup/src
walter at geodynamics.org
walter at geodynamics.org
Wed Jul 18 11:37:30 PDT 2007
Author: walter
Date: 2007-07-18 11:37:29 -0700 (Wed, 18 Jul 2007)
New Revision: 7696
Added:
long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c
long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.h
long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.meta
long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c
long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.h
long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.meta
long/3D/Gale/trunk/src/Gale/Utils/src/Wall.c
Removed:
long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c
long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.h
long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.meta
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.c
long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.h
long/3D/Gale/trunk/src/Gale/Utils/src/Init.c
long/3D/Gale/trunk/src/Gale/Utils/src/SConscript
long/3D/Gale/trunk/src/Gale/Utils/src/Utils.h
long/3D/Gale/trunk/src/Gale/Utils/src/types.h
long/3D/Gale/trunk/src/PICellerator/Utils/src/Init.c
long/3D/Gale/trunk/src/PICellerator/Utils/src/SConscript
long/3D/Gale/trunk/src/PICellerator/Utils/src/Utils.h
long/3D/Gale/trunk/src/PICellerator/Utils/src/types.h
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
Log:
r1880 at earth: boo | 2007-07-18 11:33:07 -0700
First pass at Kinetic friction. I have only tested that I did not mess up StressBC in the process
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1879
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1880
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.c 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.c 2007-07-18 18:37:29 UTC (rev 7696)
@@ -37,13 +37,7 @@
const Type FrictionVC_Type = "FrictionVC";
const Name defaultFrictionVCName = "defaultFrictionVCName";
-const char* FrictionVC_WallEnumToStr[FrictionVC_Wall_Size] = {
- "back",
- "left",
- "bottom",
- "right",
- "top",
- "front" };
+extern const char* WallEnumToStr[Wall_Size];
/*--------------------------------------------------------------------------------------------------------------------------
@@ -56,11 +50,11 @@
Dictionary* dictionary,
void* data )
{
- return (VariableCondition*)FrictionVC_New( defaultFrictionVCName, NULL, variable_Register, conFunc_Register, dictionary, (Mesh*)data );
+ return (VariableCondition*)FrictionVC_New( defaultFrictionVCName, NULL, variable_Register, conFunc_Register, dictionary, (FeMesh*)data );
}
-FrictionVC* FrictionVC_DefaultNew( Name name )
+void* FrictionVC_DefaultNew( Name name )
{
return _FrictionVC_New(
sizeof(FrictionVC),
@@ -256,8 +250,8 @@
self->isConstructed = True;
self->_dictionaryEntryName = _dictionaryEntryName;
- self->_mesh = (Mesh*)_mesh;
- self->_wall = FrictionVC_Wall_Size;
+ self->_mesh = (FeMesh*)_mesh;
+ self->_wall = Wall_Size;
self->_entryTbl = 0;
self->_entryCount = 0;
self->context=NULL;
@@ -299,20 +293,20 @@
/* Obtain which wall */
wallStr = Dictionary_Entry_Value_AsString(Dictionary_Entry_Value_GetMember(vcDictVal, "wall" ));
if (!strcasecmp(wallStr, "back"))
- self->_wall = FrictionVC_Wall_Back;
+ self->_wall = Wall_Back;
else if (!strcasecmp(wallStr, "left"))
- self->_wall = FrictionVC_Wall_Left;
+ self->_wall = Wall_Left;
else if (!strcasecmp(wallStr, "bottom"))
- self->_wall = FrictionVC_Wall_Bottom;
+ self->_wall = Wall_Bottom;
else if (!strcasecmp(wallStr, "right"))
- self->_wall = FrictionVC_Wall_Right;
+ self->_wall = Wall_Right;
else if (!strcasecmp(wallStr, "top"))
- self->_wall = FrictionVC_Wall_Top;
+ self->_wall = Wall_Top;
else if (!strcasecmp(wallStr, "front"))
- self->_wall = FrictionVC_Wall_Front;
+ self->_wall = Wall_Front;
else {
assert( 0 );
- self->_wall = FrictionVC_Wall_Size; /* invalid entry */
+ self->_wall = Wall_Size; /* invalid entry */
}
/* Obtain the variable entries */
@@ -357,7 +351,7 @@
Journal_Printf( errorStr, "Error- in %s: While parsing "
"definition of frictionVC \"%s\" (applies to wall \"%s\"), the cond. func. applied to "
"variable \"%s\" - \"%s\" - wasn't found in the c.f. register.\n",
- __func__, self->_dictionaryEntryName, FrictionVC_WallEnumToStr[self->_wall],
+ __func__, self->_dictionaryEntryName, WallEnumToStr[self->_wall],
self->_entryTbl[entry_I].varName, funcName );
Journal_Printf( errorStr, "(Available functions in the C.F. register are: ");
ConditionFunction_Register_PrintNameOfEachFunc( self->conFunc_Register, errorStr );
@@ -418,7 +412,7 @@
else
{
int i;
- self->_wall = FrictionVC_Wall_Size;
+ self->_wall = Wall_Size;
self->_entryCount = 0;
self->_entryTbl = NULL;
self->context=NULL;
@@ -464,10 +458,10 @@
Journal_Printf( info, "\t_dictionaryEntryName (ptr): %p\n", self->_dictionaryEntryName);
if (self->_dictionaryEntryName)
Journal_Printf( info, "\t\t_dictionaryEntryName: %s\n", self->_dictionaryEntryName);
- Journal_Printf( info, "\t_wall: %s\n", self->_wall == FrictionVC_Wall_Front ? "Front" :
- self->_wall == FrictionVC_Wall_Back ? "Back" : self->_wall == FrictionVC_Wall_Left ? "Left" :
- self->_wall == FrictionVC_Wall_Right ? "Right" : self->_wall == FrictionVC_Wall_Top ? "Top" :
- self->_wall == FrictionVC_Wall_Bottom ? "Bottom" : "None");
+ Journal_Printf( info, "\t_wall: %s\n", self->_wall == Wall_Front ? "Front" :
+ self->_wall == Wall_Back ? "Back" : self->_wall == Wall_Left ? "Left" :
+ self->_wall == Wall_Right ? "Right" : self->_wall == Wall_Top ? "Top" :
+ self->_wall == Wall_Bottom ? "Bottom" : "None");
Journal_Printf( info, "\t_entryCount: %u\n", self->_entryCount);
Journal_Printf( info, "\t_entryTbl (ptr): %p\n", self->_entryTbl);
if (self->_entryTbl)
@@ -557,7 +551,7 @@
newFrictionVC->friction = self->friction;
if( deep ) {
- newFrictionVC->_mesh = (Mesh*)Stg_Class_Copy( self->_mesh, NULL, deep, nameExt, map );
+ newFrictionVC->_mesh = (FeMesh*)Stg_Class_Copy( self->_mesh, NULL, deep, nameExt, map );
if( (newFrictionVC->_entryTbl = PtrMap_Find( map, self->_entryTbl )) == NULL && self->_entryTbl ) {
newFrictionVC->_entryTbl = Memory_Alloc_Array( FrictionVC_Entry, newFrictionVC->_entryCount, "FrictionVC->_entryTbl");
@@ -646,51 +640,8 @@
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;
- }
-
+ FrictionVC_get_basis_vectors(self->_wall,grid->sizes,
+ &direction,&boundary,basis);
nNodes = Mesh_GetDomainSize( self->_mesh, MT_VERTEX );
set = IndexSet_New( nNodes );
for( n_i = 0; n_i < nNodes; n_i++ ) {
@@ -699,217 +650,17 @@
( self->_mesh, MT_VERTEX, n_i ), ijk );
if(ijk[direction]==boundary)
{
- double str[6], p, normal_stress, tangential_norm;
- FeVariable_GetValueAtNode(deviatoric_stress,n_i,str);
- FeVariable_GetValueAtNode(pressure,n_i,&p);
-
- if(nDims==2)
- {
- 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)
- {
- if(self->include_lower_boundary[basis[0]]==True)
- {
- x1[0]=self->_mesh->verts[n_i][0];
- x1[1]=self->_mesh->verts[n_i][1];
- }
- else
- continue;
- }
- 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)
- {
- if(self->include_upper_boundary[basis[0]]==True)
- {
- x2[0]=self->_mesh->verts[n_i][0];
- x2[1]=self->_mesh->verts[n_i][1];
- }
- else
- continue;
- }
- 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);
+ double normal_stress, tangential_norm, n[3];
- /* Normal stress = n . stress . n */
-
- /* The ordering for components of stress is: xx,
- yy, xy */
-
- normal_stress=n[0]*n[0]*str[0]
- + n[1]*n[1]*str[1]
- + 2*n[1]*n[0]*str[2];
-
- /* 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
- {
- 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)
- {
- if(self->include_lower_boundary[basis[0]]==True)
- {
- x1[0]=self->_mesh->verts[n_i][0];
- x1[1]=self->_mesh->verts[n_i][1];
- }
- else
- continue;
- }
- 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)
- {
- if(self->include_upper_boundary[basis[0]]==True)
- {
- x2[0]=self->_mesh->verts[n_i][0];
- x2[1]=self->_mesh->verts[n_i][1];
- }
- else
- continue;
- }
- 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)
- {
- if(self->include_lower_boundary[basis[1]]==True)
- {
- x3[0]=self->_mesh->verts[n_i][0];
- x3[1]=self->_mesh->verts[n_i][1];
- }
- else
- continue;
- }
- 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)
- {
- if(self->include_upper_boundary[basis[1]]==True)
- {
- x4[0]=self->_mesh->verts[n_i][0];
- x4[1]=self->_mesh->verts[n_i][1];
- }
- else
- continue;
- }
- 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);
- }
-
/* Finally, add in the point if the frictional force
keeps the material pinned to the boundary. */
-/* printf("Adding %d %d %lf %lf %lf\n",ijk[0],ijk[1],p,normal_stress,tangential_norm); */
- if(self->friction*(p-normal_stress)>=tangential_norm)
+ if(FrictionVC_get_interface_stresses
+ (pressure,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*normal_stress>=tangential_norm)
{
IndexSet_Add(set,n_i);
}
@@ -919,7 +670,7 @@
else
{
switch (self->_wall) {
- case FrictionVC_Wall_Front:
+ case Wall_Front:
if ( nDims < 3 || !grid->sizes[2] ) {
set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
@@ -928,7 +679,7 @@
}
break;
- case FrictionVC_Wall_Back:
+ case Wall_Back:
if ( nDims < 3 || !grid->sizes[2] ) {
set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
@@ -937,7 +688,7 @@
}
break;
- case FrictionVC_Wall_Top:
+ case Wall_Top:
if ( nDims < 2 || !grid->sizes[1] ) {
set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
@@ -946,7 +697,7 @@
}
break;
- case FrictionVC_Wall_Bottom:
+ case Wall_Bottom:
if ( nDims < 2 || !grid->sizes[1] ) {
set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
@@ -955,7 +706,7 @@
}
break;
- case FrictionVC_Wall_Left:
+ case Wall_Left:
if ( nDims < 1 ) {
set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
@@ -964,7 +715,7 @@
}
break;
- case FrictionVC_Wall_Right:
+ case Wall_Right:
if( nDims < 1 ) {
set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
@@ -973,7 +724,7 @@
}
break;
- case FrictionVC_Wall_Size:
+ case Wall_Size:
default:
assert(0);
break;
@@ -1038,12 +789,12 @@
Journal_Printf( stream, "\ttype: %s, set: ", self->type );
Journal_Printf( stream, "%s\n",
- self->_wall == FrictionVC_Wall_Front ? "Front" :
- self->_wall == FrictionVC_Wall_Back ? "Back" :
- self->_wall == FrictionVC_Wall_Left ? "Left" :
- self->_wall == FrictionVC_Wall_Right ? "Right" :
- self->_wall == FrictionVC_Wall_Top ? "Top" :
- self->_wall == FrictionVC_Wall_Bottom ? "Bottom" : "None" );
+ self->_wall == Wall_Front ? "Front" :
+ self->_wall == Wall_Back ? "Back" :
+ self->_wall == Wall_Left ? "Left" :
+ self->_wall == Wall_Right ? "Right" :
+ self->_wall == Wall_Top ? "Top" :
+ self->_wall == Wall_Bottom ? "Bottom" : "None" );
}
/*--------------------------------------------------------------------------------------------------------------------------
@@ -1054,3 +805,275 @@
/*--------------------------------------------------------------------------------------------------------------------------
** Functions
*/
+
+void FrictionVC_get_basis_vectors(Wall wall, unsigned int sizes[],
+ unsigned *direction, unsigned *boundary,
+ unsigned basis[])
+{
+ /* Note that we select the order of the basis vectors to
+ point inwards */
+ switch (wall) {
+ case Wall_Front:
+ *direction=K_AXIS;
+ *boundary=sizes[*direction]-1;
+ basis[0]=J_AXIS;
+ basis[1]=I_AXIS;
+ break;
+ case Wall_Back:
+ *direction=K_AXIS;
+ *boundary=0;
+ basis[0]=I_AXIS;
+ basis[1]=J_AXIS;
+ break;
+ case Wall_Top:
+ *direction=J_AXIS;
+ *boundary=sizes[*direction]-1;
+ basis[0]=I_AXIS;
+ basis[1]=K_AXIS;
+ break;
+ case Wall_Bottom:
+ *direction=J_AXIS;
+ *boundary=0;
+ basis[0]=K_AXIS;
+ basis[1]=I_AXIS;
+ break;
+ case Wall_Left:
+ *direction=I_AXIS;
+ *boundary=0;
+ basis[0]=K_AXIS;
+ basis[1]=J_AXIS;
+ break;
+ case Wall_Right:
+ *direction=I_AXIS;
+ *boundary=sizes[*direction]-1;
+ basis[0]=J_AXIS;
+ basis[1]=K_AXIS;
+ break;
+ case Wall_Size:
+ default:
+ assert(0);
+ break;
+ }
+}
+
+/* This function returns the normal to the surface and the magnitude
+ of the normal and tangential stresses. The exit status tells if
+ the point is on a non-included boundary. */
+
+Bool FrictionVC_get_interface_stresses(FeVariable *pressure,
+ FeVariable *deviatoric_stress,
+ FeMesh *mesh,
+ unsigned sizes[],
+ Node_LocalIndex n_i,
+ IJK ijk,
+ unsigned nDims,
+ unsigned basis[],
+ Bool include_lower_boundary[],
+ Bool include_upper_boundary[],
+ double *normal_stress,
+ double *tangential_norm,
+ double n[])
+{
+ double str[6], p;
+ FeVariable_GetValueAtNode(deviatoric_stress,n_i,str);
+ FeVariable_GetValueAtNode(pressure,n_i,&p);
+
+ if(nDims==2)
+ {
+ double x1[2],x2[2],surface[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)
+ {
+ if(include_lower_boundary[basis[0]]==True)
+ {
+ x1[0]=mesh->verts[n_i][0];
+ x1[1]=mesh->verts[n_i][1];
+ }
+ else
+ return False;
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[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(include_upper_boundary[basis[0]]==True)
+ {
+ x2[0]=mesh->verts[n_i][0];
+ x2[1]=mesh->verts[n_i][1];
+ }
+ else
+ return False;
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[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];
+ }
+ /* Get the surface */
+ Vec_Sub2D(surface,x2,x1);
+ n[0]=sign*surface[1];
+ n[1]=sign*surface[0];
+ Vec_Norm2D(n,n);
+
+ /* Normal stress = n . stress . n */
+
+ /* The ordering for components of stress is: xx,
+ yy, xy */
+
+ *normal_stress=n[0]*n[0]*str[0]
+ + n[1]*n[1]*str[1]
+ + 2*n[1]*n[0]*str[2];
+
+ /* 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
+ {
+ double x1[3],x2[3],x3[3],x4[3],line1[3],line2[3],tangential_stress[3];
+ unsigned n_temp, ijk_temp[3];
+
+ /* Get x1 */
+ if(ijk[basis[0]]==0)
+ {
+ if(include_lower_boundary[basis[0]]==True)
+ {
+ x1[0]=mesh->verts[n_i][0];
+ x1[1]=mesh->verts[n_i][1];
+ }
+ else
+ return False;
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[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(include_upper_boundary[basis[0]]==True)
+ {
+ x2[0]=mesh->verts[n_i][0];
+ x2[1]=mesh->verts[n_i][1];
+ }
+ else
+ return False;
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[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];
+ }
+
+ /* Get x3 */
+ if(ijk[basis[1]]==0)
+ {
+ if(include_lower_boundary[basis[1]]==True)
+ {
+ x3[0]=mesh->verts[n_i][0];
+ x3[1]=mesh->verts[n_i][1];
+ }
+ else
+ return False;
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[basis[1]]-=1;
+ n_temp=RegularMeshUtils_Node_3DTo1D(mesh,ijk_temp);
+ x3[0]=mesh->verts[n_temp][0];
+ x3[1]=mesh->verts[n_temp][1];
+ }
+ /* Get x4 */
+ if(ijk[basis[1]]==sizes[basis[1]]-1)
+ {
+ if(include_upper_boundary[basis[1]]==True)
+ {
+ x4[0]=mesh->verts[n_i][0];
+ x4[1]=mesh->verts[n_i][1];
+ }
+ else
+ return False;
+ }
+ else
+ {
+ Vec_Set2D(ijk_temp,ijk);
+ ijk_temp[basis[1]]+=1;
+ n_temp=RegularMeshUtils_Node_3DTo1D(mesh,ijk_temp);
+ x4[0]=mesh->verts[n_temp][0];
+ x4[1]=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);
+ }
+ /* Add in the pressure. The sign convention switches the direction
+ of stress from the deviatoric stress. */
+ *normal_stress=p-(*normal_stress);
+}
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.h 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/FrictionVC.h 2007-07-18 18:37:29 UTC (rev 7696)
@@ -36,7 +36,7 @@
extern const Type FrictionVC_Type;
- extern const char* FrictionVC_WallEnumToStr[FrictionVC_Wall_Size];
+ extern const char* WallEnumToStr[Wall_Size];
#define __FrictionVC_Entry \
Name varName; \
@@ -53,10 +53,10 @@
\
/* Stg_Class info */ \
Name _dictionaryEntryName; \
- FrictionVC_Wall _wall; \
+ Wall _wall; \
FrictionVC_Entry_Index _entryCount; \
FrictionVC_Entry* _entryTbl; \
- Mesh* _mesh; \
+ FeMesh* _mesh; \
/* I would like to make this a FiniteElementContext*, but */ \
/* then there are problems compiling because this is in */ \
/* StGermain, and we do not have access to StgFEM yet. */ \
@@ -79,7 +79,7 @@
Dictionary* dictionary,
void* data );
- FrictionVC* FrictionVC_DefaultNew( Name name );
+ void* FrictionVC_DefaultNew( Name name );
FrictionVC* FrictionVC_New(
Name name,
@@ -196,5 +196,22 @@
** Functions
*/
+Bool FrictionVC_get_interface_stresses(FeVariable *pressure,
+ FeVariable *deviatoric_stress,
+ FeMesh *mesh,
+ unsigned sizes[],
+ Node_LocalIndex n_i,
+ IJK ijk,
+ unsigned nDims,
+ unsigned basis[],
+ Bool include_lower_boundary[],
+ Bool include_upper_boundary[],
+ double *normal_stress,
+ double *tangential_norm,
+ double n[]);
+
+void FrictionVC_get_basis_vectors(Wall wall, unsigned int sizes[],
+ unsigned *direction, unsigned *boundary,
+ unsigned basis[]);
#endif /* __Discretisation_Utils_FrictionVC_h__ */
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/Init.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/Init.c 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/Init.c 2007-07-18 18:37:29 UTC (rev 7696)
@@ -57,13 +57,22 @@
Journal_Printf( Journal_Register( DebugStream_Type, "Context" ), "In: %s\n", __func__ ); /* DO NOT CHANGE OR REMOVE */
- Stg_ComponentRegister_Add( componentRegister, GaleContext_Type, "0", GaleContext_DefaultNew );
+ Stg_ComponentRegister_Add
+ ( componentRegister, GaleContext_Type, "0", GaleContext_DefaultNew );
+
RegisterParent( GaleContext_Type, UnderworldContext_Type );
VariableCondition_Register_Add( variableCondition_Register, FrictionVC_Type, FrictionVC_Factory );
Stg_ComponentRegister_Add( Stg_ComponentRegister_Get_ComponentRegister(), FrictionVC_Type,
- "0", (void* (*)(Name))FrictionVC_DefaultNew );
+ "0", FrictionVC_DefaultNew );
RegisterParent( FrictionVC_Type, VariableCondition_Type );
+ Stg_ComponentRegister_Add( componentRegister, StressBC_Type, "0", _StressBC_DefaultNew );
+ RegisterParent( StressBC_Type, ForceTerm_Type );
+ Stg_ComponentRegister_Add( componentRegister, KineticFriction_Type,
+ "0", _KineticFriction_DefaultNew );
+
+ RegisterParent( KineticFriction_Type, ForceTerm_Type );
+
return True;
}
Added: long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.c 2007-07-18 18:37:29 UTC (rev 7696)
@@ -0,0 +1,464 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+** Melbourne, 3053, Australia.
+** Copyright (c) 2005-2006, Monash Cluster Computing, Building 28, Monash University Clayton Campus,
+** Victoria, 3800, Australia
+**
+** Primary Contributing Organisations:
+** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+** Australian Computational Earth Systems Simulator - http://www.access.edu.au
+** Monash Cluster Computing - http://www.mcc.monash.edu.au
+**
+** Contributors:
+** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+** Patrick D. Sunter, Software Engineer, VPAC. (patrick at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** David Stegman, Postdoctoral Fellow, Monash University. (david.stegman at sci.monash.edu.au)
+** Wendy Sharples, PhD Student, Monash University (wendy.sharples at sci.monash.edu.au)
+**
+** Copyright (C) 2006, 2007 California Institute of Technology
+** by Walter Landry
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** $Id$
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StgFEM/StgFEM.h>
+#include <PICellerator/Voronoi/Voronoi.h>
+#include <PICellerator/PopulationControl/PopulationControl.h>
+#include <PICellerator/Weights/Weights.h>
+#include <PICellerator/MaterialPoints/MaterialPoints.h>
+
+#include "PICellerator/Utils/types.h"
+#include "PICellerator/Utils/MaterialSwarmVariable.h"
+
+#include "types.h"
+#include "KineticFriction.h"
+#include "FrictionVC.h"
+
+#include <assert.h>
+#include <string.h>
+
+/* Textual name of this class */
+const Type KineticFriction_Type = "KineticFriction";
+
+extern const char* WallEnumToStr[];
+
+KineticFriction* KineticFriction_New
+(
+ Name name,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ KineticFriction* self = (KineticFriction*) _KineticFriction_DefaultNew( name );
+
+ KineticFriction_InitAll
+ (
+ self,
+ forceVector,
+ integrationSwarm,
+ conFunc_Register,
+ context);
+
+ return self;
+}
+
+/* Creation implementation / Virtual constructor */
+KineticFriction* _KineticFriction_New
+(
+ SizeT sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ ForceTerm_AssembleElementFunction* _assembleElement,
+ ConditionFunction_Register* conFunc_Register,
+ Name name,
+ FiniteElementContext* context)
+{
+ KineticFriction* self;
+
+ /* Allocate memory */
+ assert( sizeOfSelf >= sizeof(KineticFriction) );
+ self = (KineticFriction*) _ForceTerm_New
+ (
+ sizeOfSelf,
+ type,
+ _delete,
+ _print,
+ _copy,
+ _defaultConstructor,
+ _construct,
+ _build,
+ _initialise,
+ _execute,
+ _destroy,
+ _assembleElement,
+ name );
+
+ self->conFunc_Register=conFunc_Register;
+ self->context=context;
+
+ return self;
+}
+
+void _KineticFriction_Init
+(
+ KineticFriction* self,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ self->conFunc_Register=conFunc_Register;
+ self->context=context;
+ self->pressure_name=NULL;
+ self->deviatoric_stress_name=NULL;
+}
+
+void KineticFriction_InitAll
+(
+ void* forceTerm,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ KineticFriction* self = (KineticFriction*) forceTerm;
+
+ ForceTerm_InitAll( self, forceVector, integrationSwarm, NULL );
+ _KineticFriction_Init( self, conFunc_Register, context );
+}
+
+void _KineticFriction_Delete( void* forceTerm ) {
+ KineticFriction* self = (KineticFriction*)forceTerm;
+
+ _ForceTerm_Delete( self );
+}
+
+void _KineticFriction_Print( void* forceTerm, Stream* stream ) {
+ KineticFriction* self = (KineticFriction*)forceTerm;
+
+ _ForceTerm_Print( self, stream );
+}
+
+void* _KineticFriction_DefaultNew( Name name ) {
+ return (void*)_KineticFriction_New
+ (
+ sizeof(KineticFriction),
+ KineticFriction_Type,
+ _KineticFriction_Delete,
+ _KineticFriction_Print,
+ NULL,
+ _KineticFriction_DefaultNew,
+ _KineticFriction_Construct,
+ _KineticFriction_Build,
+ _KineticFriction_Initialise,
+ _KineticFriction_Execute,
+ _KineticFriction_Destroy,
+ _KineticFriction_AssembleElement,
+ NULL,
+ name,
+ NULL);
+}
+
+void _KineticFriction_Construct( void* forceTerm, Stg_ComponentFactory* cf,
+ void* data ) {
+ KineticFriction* self = (KineticFriction*)forceTerm;
+ FiniteElementContext* context;
+
+ /* Construct Parent */
+ _ForceTerm_Construct( self, cf, data );
+
+ context = (FiniteElementContext*)Stg_ComponentFactory_ConstructByName
+ ( cf, "context", FiniteElementContext, True, data ) ;
+
+ _KineticFriction_Init( self, context->condFunc_Register, context );
+ {
+ char* wallStr;
+
+ /* Obtain which wall */
+ wallStr = Stg_ComponentFactory_GetString( cf, self->name,
+ "wall", "");
+ if (!strcasecmp(wallStr, "back"))
+ self->_wall = Wall_Back;
+ else if (!strcasecmp(wallStr, "left"))
+ self->_wall = Wall_Left;
+ else if (!strcasecmp(wallStr, "bottom"))
+ self->_wall = Wall_Bottom;
+ else if (!strcasecmp(wallStr, "right"))
+ self->_wall = Wall_Right;
+ else if (!strcasecmp(wallStr, "top"))
+ self->_wall = Wall_Top;
+ else if (!strcasecmp(wallStr, "front"))
+ self->_wall = Wall_Front;
+ else {
+ assert( 0 );
+ self->_wall = Wall_Size; /* invalid entry */
+ }
+ }
+
+ self->pressure_name = Stg_ComponentFactory_GetString( cf, self->name,
+ "PressureField", "");
+ self->deviatoric_stress_name =
+ Stg_ComponentFactory_GetString( cf, self->name, "StressField", "");
+ self->velocity_name = Stg_ComponentFactory_GetString( cf, self->name,
+ "VelocityField", "");
+
+ _KineticFriction_GetValues(cf,self,&(self->v_entry[I_AXIS]),"vx");
+ _KineticFriction_GetValues(cf,self,&(self->v_entry[J_AXIS]),"vy");
+ _KineticFriction_GetValues(cf,self,&(self->v_entry[K_AXIS]),"vz");
+ _KineticFriction_GetValues(cf,self,&(self->friction_entry),"friction");
+
+}
+
+/* Gets the actual values used by KineticFriction (e.g. a float or
+ a function). */
+void _KineticFriction_GetValues(Stg_ComponentFactory* cf,
+ void *kineticFriction,
+ StressBC_Entry *entry, char *prefix)
+{
+ KineticFriction* self = (KineticFriction*)kineticFriction;
+ char temp_str[20];
+ char *type;
+
+ strcat(strcpy(temp_str,prefix),"_type");
+ type=Stg_ComponentFactory_GetString( cf, self->name, temp_str, "");
+ strcat(strcpy(temp_str,prefix),"_value");
+
+ if(!strcasecmp(type,"double") || !strcasecmp(type,"float")
+ || strlen(type)==0)
+ {
+ entry->type = StressBC_Double;
+ entry->DoubleValue =
+ Stg_ComponentFactory_GetDouble( cf, self->name, temp_str, 0.0);
+ }
+ else if(!strcasecmp(type,"func"))
+ {
+ char *funcName =
+ Stg_ComponentFactory_GetString( cf, self->name, temp_str, "");
+ Index cfIndex;
+
+ cfIndex=ConditionFunction_Register_GetIndex(self->conFunc_Register,
+ funcName);
+ entry->type = StressBC_ConditionFunction;
+ if ( cfIndex == (unsigned)-1 ) {
+ Stream* errorStr = Journal_Register( Error_Type, self->type );
+
+ Journal_Printf( errorStr, "Error- in %s: While parsing "
+ "definition of KineticFriction (applies to wall \"%s\"), the cond. func. applied to "
+ "\"%s\" - \"%s\" - wasn't found in the c.f. register.\n",
+ __func__, WallEnumToStr[self->_wall],
+ prefix, funcName );
+ Journal_Printf( errorStr, "(Available functions in the C.F. register are: ");
+ ConditionFunction_Register_PrintNameOfEachFunc
+ ( self->conFunc_Register, errorStr );
+ Journal_Printf( errorStr, ")\n");
+ assert(0);
+ }
+ entry->CFIndex = cfIndex;
+ }
+ else if(strlen(type)!=0)
+ {
+ Stream* errorStr = Journal_Register( Error_Type, self->type );
+ Journal_Printf( errorStr, "Error- in %s: While parsing "
+ "definition of KineticFriction (applies to wall \"%s\"), the type of condition \"%s\"\nis not supported. Supported types are \"double\" and \"function\".",
+ __func__, WallEnumToStr[self->_wall],
+ type );
+ assert(0);
+ }
+}
+
+void _KineticFriction_Build( void* forceTerm, void* data ) {
+ KineticFriction* self = (KineticFriction*)forceTerm;
+ Name name;
+
+ _ForceTerm_Build( self, data );
+}
+
+void _KineticFriction_Initialise( void* forceTerm, void* data ) {
+ KineticFriction* self = (KineticFriction*)forceTerm;
+ Index i;
+
+ _ForceTerm_Initialise( self, data );
+
+}
+
+void _KineticFriction_Execute( void* forceTerm, void* data ) {
+ _ForceTerm_Execute( forceTerm, data );
+}
+
+void _KineticFriction_Destroy( void* forceTerm, void* data ) {
+ _ForceTerm_Destroy( forceTerm, data );
+}
+
+void _KineticFriction_AssembleElement( void* forceTerm,
+ ForceVector* forceVector,
+ Element_LocalIndex lElement_I,
+ double* elForceVec ) {
+ KineticFriction* self = (KineticFriction*) forceTerm;
+ Element_NodeIndex elementNodeCount;
+ Dimension_Index dim = forceVector->dim;
+ FeMesh* mesh = forceVector->feVariable->feMesh;
+ Grid* grid;
+ Node_ElementLocalIndex eNode_I;
+ ElementType* elementType;
+ Dof_Index nodeDofCount;
+ double stress, area, n[3], v_diff[3], v_boundary[3], v_dot_n,
+ v_perpendicular[3], v_norm, friction_stress;
+ IJK ijk;
+ unsigned direction, boundary, basis[2], d;
+ FeVariable *pressure, *deviatoric_stress, *velocity;
+ pressure=deviatoric_stress=velocity=NULL;
+
+ Node_DomainIndex* elementNodes = NULL;
+
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ nodeDofCount = dim;
+
+ grid = *(Grid**)ExtensionManager_Get
+ ( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+
+ area=StressBC_compute_face_area(self->_wall,mesh,lElement_I,
+ &elementNodeCount,dim,&elementNodes);
+
+ /* Get the basis vectors for the boundary */
+ FrictionVC_get_basis_vectors(self->_wall,grid->sizes,
+ &direction,&boundary,basis);
+
+ /* Apply the stress */
+ for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
+ /* Make sure that we are on the boundary */
+ int entry_I;
+ ConditionFunction* cf;
+ RegularMeshUtils_Node_1DTo3D
+ (mesh,Mesh_DomainToGlobal(mesh,MT_VERTEX,elementNodes[eNode_I]),ijk);
+
+ if(ijk[direction]==boundary)
+ {
+ unsigned overcount=StressBC_get_overcount(ijk,grid->sizes);
+ double v[3], v_surface[3], v_surface_diff[3], normal_stress,
+ tangential_norm, friction_coefficient;
+ Bool include_boundary[3];
+ 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[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+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 StressnVC, 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
+ FrictionVC's. */
+ for(d=0;d<dim;d++)
+ include_boundary[d]=True;
+ FrictionVC_get_interface_stresses
+ (pressure,deviatoric_stress,mesh,grid->sizes,elementNodes[eNode_I],
+ ijk,dim,basis,include_boundary,
+ include_boundary,&normal_stress,&tangential_norm,n);
+
+ 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,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);
+ }
+ friction_stress=friction_coefficient*normal_stress*area/overcount;
+ for(d=0;d<dim;++d)
+ elForceVec[ eNode_I*nodeDofCount + d]+=v_norm*friction_stress;
+ }
+ }
+}
+
Added: long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.h 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.h 2007-07-18 18:37:29 UTC (rev 7696)
@@ -0,0 +1,115 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+** Melbourne, 3053, Australia.
+** Copyright (c) 2005-2006, Monash Cluster Computing, Building 28, Monash University Clayton Campus,
+** Victoria, 3800, Australia
+**
+** Primary Contributing Organisations:
+** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+** Australian Computational Earth Systems Simulator - http://www.access.edu.au
+** Monash Cluster Computing - http://www.mcc.monash.edu.au
+**
+** Contributors:
+** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+** Patrick D. Sunter, Software Engineer, VPAC. (patrick at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** David Stegman, Postdoctoral Fellow, Monash University. (david.stegman at sci.monash.edu.au)
+** Wendy Sharples, PhD Student, Monash University (wendy.sharples at sci.monash.edu.au)
+**
+** Modified 2006 by Walter Landry (California Institute of Technology)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+*/
+
+
+#ifndef __Gale_Utils_KineticFriction_h__
+#define __Gale_Utils_KineticFriction_h__
+
+ /** Textual name of this class */
+ extern const Type KineticFriction_Type;
+
+ /** KineticFriction class contents */
+ #define __KineticFriction \
+ /* General info */ \
+ __ForceTerm \
+ \
+ /* KineticFriction info */ \
+ Wall _wall; \
+ StressBC_Entry v_entry[3]; \
+ StressBC_Entry friction_entry; \
+ ConditionFunction_Register* conFunc_Register; \
+ FiniteElementContext* context; \
+ char *deviatoric_stress_name; \
+ char *pressure_name; \
+ char *velocity_name; \
+
+ struct KineticFriction { __KineticFriction };
+
+ KineticFriction* KineticFriction_New(
+ Name name,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context);
+
+ KineticFriction* _KineticFriction_New(
+ SizeT sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ ForceTerm_AssembleElementFunction* _assembleElement,
+ ConditionFunction_Register* conFunc_Register,
+ Name name,
+ FiniteElementContext* context);
+
+ void KineticFriction_InitAll(
+ void* forceTerm,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context);
+
+ void _KineticFriction_Delete( void* forceTerm );
+ void _KineticFriction_Print( void* forceTerm, Stream* stream );
+
+ void* _KineticFriction_DefaultNew( Name name ) ;
+ void _KineticFriction_Construct( void* forceTerm, Stg_ComponentFactory* cf, void* data ) ;
+ void _KineticFriction_Build( void* forceTerm, void* data ) ;
+ void _KineticFriction_Initialise( void* forceTerm, void* data ) ;
+ void _KineticFriction_Execute( void* forceTerm, void* data ) ;
+ void _KineticFriction_Destroy( void* forceTerm, void* data ) ;
+
+ void _KineticFriction_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) ;
+ void _KineticFriction_GetValues(Stg_ComponentFactory* cf,
+ void *kineticFriction,
+ StressBC_Entry *entry,
+ char *prefix);
+#endif
Added: long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.meta
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.meta 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/KineticFriction.meta 2007-07-18 18:37:29 UTC (rev 7696)
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">KineticFriction</param>
+<param name="Organisation">VPAC and MCC</param>
+<param name="Project">Gale</param>
+<param name="Location">./Gale/Utils/src/</param>
+<param name="Project Web">https://csd.vpac.org/twiki/bin/view/PICellerator/WebHome</param>
+<param name="Copyright">Copyright (C) 2005 VPAC and Monash Cluster Computing.</param>
+<param name="License">https://csd.vpac.org/twiki/bin/view/Stgermain/SoftwareLicense http://www.opensource.org/licenses/bsd-license.php</param>
+<param name="Parent">ForceTerm</param>
+<param name="Description">This adds a stress boundary condition.</param>
+
+<!--Now the interesting stuff-->
+
+
+<list name="Params">
+ <struct>
+ </struct>
+
+</list>
+
+<!-- Add an exmaple XML if possible -->
+<param name="Example">
+ <struct name="stressBC">
+ <param name="Type">KineticFriction</param>
+ <param name="ForceVector">mom_force</param>
+ <param name="Swarm">materialPoints</param>
+ <param name="wall">bottom</param>
+ <param name="friction_type">double</param>
+ <param name="friction_value">1.0</param>
+ </struct>
+</param>
+
+</StGermainData>
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/SConscript
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/SConscript 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/SConscript 2007-07-18 18:37:29 UTC (rev 7696)
@@ -23,6 +23,8 @@
Finalise.h
FrictionVC.h
Init.h
+KineticFriction.h
+StressBC.h
types.h
Utils.h""")
@@ -30,10 +32,17 @@
c_files=Split("""GaleContext.c
Finalise.c
FrictionVC.c
-Init.c""")
+Init.c
+KineticFriction.c
+StressBC.c
+Wall.c
+""")
-meta_files=Split("""GaleContext.meta FrictionVC.meta""")
+meta_files=Split("""GaleContext.meta
+FrictionVC.meta
+KineticFriction.meta
+StressBC.meta""")
c_files+=[local_env.Meta(meta_files)]
Added: long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c 2007-07-18 18:37:29 UTC (rev 7696)
@@ -0,0 +1,531 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+** Melbourne, 3053, Australia.
+** Copyright (c) 2005-2006, Monash Cluster Computing, Building 28, Monash University Clayton Campus,
+** Victoria, 3800, Australia
+**
+** Primary Contributing Organisations:
+** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+** Australian Computational Earth Systems Simulator - http://www.access.edu.au
+** Monash Cluster Computing - http://www.mcc.monash.edu.au
+**
+** Contributors:
+** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+** Patrick D. Sunter, Software Engineer, VPAC. (patrick at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** David Stegman, Postdoctoral Fellow, Monash University. (david.stegman at sci.monash.edu.au)
+** Wendy Sharples, PhD Student, Monash University (wendy.sharples at sci.monash.edu.au)
+**
+** Modified 2006 by Walter Landry (California Institute of Technology)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** $Id$
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StgFEM/StgFEM.h>
+#include <PICellerator/Voronoi/Voronoi.h>
+#include <PICellerator/PopulationControl/PopulationControl.h>
+#include <PICellerator/Weights/Weights.h>
+#include <PICellerator/MaterialPoints/MaterialPoints.h>
+
+#include "PICellerator/Utils/types.h"
+#include "PICellerator/Utils/MaterialSwarmVariable.h"
+
+#include "types.h"
+#include "StressBC.h"
+
+#include <assert.h>
+#include <string.h>
+
+/* Textual name of this class */
+const Type StressBC_Type = "StressBC";
+
+extern const char* WallEnumToStr[Wall_Size];
+
+StressBC* StressBC_New(
+ Name name,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ StressBC* self = (StressBC*) _StressBC_DefaultNew( name );
+
+ StressBC_InitAll(
+ self,
+ forceVector,
+ integrationSwarm,
+ conFunc_Register,
+ context);
+
+ return self;
+}
+
+/* Creation implementation / Virtual constructor */
+StressBC* _StressBC_New(
+ SizeT sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ ForceTerm_AssembleElementFunction* _assembleElement,
+ ConditionFunction_Register* conFunc_Register,
+ Name name,
+ FiniteElementContext* context)
+{
+ StressBC* self;
+
+ /* Allocate memory */
+ assert( sizeOfSelf >= sizeof(StressBC) );
+ self = (StressBC*) _ForceTerm_New(
+ sizeOfSelf,
+ type,
+ _delete,
+ _print,
+ _copy,
+ _defaultConstructor,
+ _construct,
+ _build,
+ _initialise,
+ _execute,
+ _destroy,
+ _assembleElement,
+ name );
+
+ self->conFunc_Register=conFunc_Register;
+ self->context=context;
+
+ return self;
+}
+
+void _StressBC_Init(
+ StressBC* self,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ self->numEntries = 0;
+ self->conFunc_Register=conFunc_Register;
+ self->context=context;
+}
+
+void StressBC_InitAll(
+ void* forceTerm,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ StressBC* self = (StressBC*) forceTerm;
+
+ ForceTerm_InitAll( self, forceVector, integrationSwarm, NULL );
+ _StressBC_Init( self, conFunc_Register, context );
+}
+
+void _StressBC_Delete( void* forceTerm ) {
+ StressBC* self = (StressBC*)forceTerm;
+
+ _ForceTerm_Delete( self );
+}
+
+void _StressBC_Print( void* forceTerm, Stream* stream ) {
+ StressBC* self = (StressBC*)forceTerm;
+
+ _ForceTerm_Print( self, stream );
+}
+
+void* _StressBC_DefaultNew( Name name ) {
+ return (void*)_StressBC_New(
+ sizeof(StressBC),
+ StressBC_Type,
+ _StressBC_Delete,
+ _StressBC_Print,
+ NULL,
+ _StressBC_DefaultNew,
+ _StressBC_Construct,
+ _StressBC_Build,
+ _StressBC_Initialise,
+ _StressBC_Execute,
+ _StressBC_Destroy,
+ _StressBC_AssembleElement,
+ NULL,
+ name,
+ NULL);
+}
+
+void _StressBC_Construct( void* forceTerm, Stg_ComponentFactory* cf, void* data ) {
+ StressBC* self = (StressBC*)forceTerm;
+ FiniteElementContext* context;
+
+ /* Construct Parent */
+ _ForceTerm_Construct( self, cf, data );
+
+ context = (FiniteElementContext*)Stg_ComponentFactory_ConstructByName
+ ( cf, "context", FiniteElementContext, True, data ) ;
+
+ _StressBC_Init( self, context->condFunc_Register, context );
+ {
+ char* wallStr;
+
+ /* Obtain which wall */
+ wallStr = Stg_ComponentFactory_GetString( cf, self->name,
+ "wall", "");
+ if (!strcasecmp(wallStr, "back"))
+ self->_wall = Wall_Back;
+ else if (!strcasecmp(wallStr, "left"))
+ self->_wall = Wall_Left;
+ else if (!strcasecmp(wallStr, "bottom"))
+ self->_wall = Wall_Bottom;
+ else if (!strcasecmp(wallStr, "right"))
+ self->_wall = Wall_Right;
+ else if (!strcasecmp(wallStr, "top"))
+ self->_wall = Wall_Top;
+ else if (!strcasecmp(wallStr, "front"))
+ self->_wall = Wall_Front;
+ else {
+ assert( 0 );
+ self->_wall = Wall_Size; /* invalid entry */
+ }
+ }
+ _StressBC_GetValues(cf,self,"x");
+ _StressBC_GetValues(cf,self,"y");
+ _StressBC_GetValues(cf,self,"z");
+
+}
+
+/* Gets the actual values used by the StressBC (e.g. a float or a function). */
+void _StressBC_GetValues(Stg_ComponentFactory* cf, void *stressBC,
+ char *direction)
+{
+ StressBC* self = (StressBC*)stressBC;
+ char temp_str[20];
+ char *type;
+
+ switch(*direction)
+ {
+ case 'x':
+ self->_entryTbl[self->numEntries].axis=I_AXIS;
+ break;
+ case 'y':
+ self->_entryTbl[self->numEntries].axis=J_AXIS;
+ break;
+ case 'z':
+ self->_entryTbl[self->numEntries].axis=K_AXIS;
+ break;
+ }
+
+ strcat(strcpy(temp_str,direction),"_type");
+ type=Stg_ComponentFactory_GetString( cf, self->name, temp_str, "");
+ strcat(strcpy(temp_str,direction),"_value");
+
+ if(!strcasecmp(type,"double") || !strcasecmp(type,"float"))
+ {
+ self->_entryTbl[self->numEntries].type = StressBC_Double;
+ self->_entryTbl[self->numEntries].DoubleValue =
+ Stg_ComponentFactory_GetDouble( cf, self->name, temp_str, 0.0);
+ (self->numEntries)++;
+ }
+ else if(!strcasecmp(type,"func"))
+ {
+ char *funcName =
+ Stg_ComponentFactory_GetString( cf, self->name, temp_str, "");
+ Index cfIndex;
+
+ cfIndex = ConditionFunction_Register_GetIndex
+ ( self->conFunc_Register, funcName);
+ self->_entryTbl[self->numEntries].type =
+ StressBC_ConditionFunction;
+ if ( cfIndex == (unsigned)-1 ) {
+ Stream* errorStr = Journal_Register( Error_Type, self->type );
+
+ Journal_Printf( errorStr, "Error- in %s: While parsing "
+ "definition of StressBC (applies to wall \"%s\"), the cond. func. applied to "
+ "direction \"%s\" - \"%s\" - wasn't found in the c.f. register.\n",
+ __func__, WallEnumToStr[self->_wall],
+ direction, funcName );
+ Journal_Printf( errorStr, "(Available functions in the C.F. register are: ");
+ ConditionFunction_Register_PrintNameOfEachFunc
+ ( self->conFunc_Register, errorStr );
+ Journal_Printf( errorStr, ")\n");
+ assert(0);
+ }
+ self->_entryTbl[self->numEntries].CFIndex = cfIndex;
+ (self->numEntries)++;
+ }
+ else if(strlen(type)!=0)
+ {
+ Stream* errorStr = Journal_Register( Error_Type, self->type );
+ Journal_Printf( errorStr, "Error- in %s: While parsing "
+ "definition of StressBC (applies to wall \"%s\"), the type of condition \"%s\"\nis not supported. Supported types are \"double\" and \"function\".",
+ __func__, WallEnumToStr[self->_wall],
+ type );
+ assert(0);
+ }
+}
+
+void _StressBC_Build( void* forceTerm, void* data ) {
+ StressBC* self = (StressBC*)forceTerm;
+ Name name;
+
+ _ForceTerm_Build( self, data );
+}
+
+void _StressBC_Initialise( void* forceTerm, void* data ) {
+ StressBC* self = (StressBC*)forceTerm;
+ Index i;
+
+ _ForceTerm_Initialise( self, data );
+
+}
+
+void _StressBC_Execute( void* forceTerm, void* data ) {
+ _ForceTerm_Execute( forceTerm, data );
+}
+
+void _StressBC_Destroy( void* forceTerm, void* data ) {
+ _ForceTerm_Destroy( forceTerm, data );
+}
+
+void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
+ StressBC* self = (StressBC*) forceTerm;
+ Element_NodeIndex elementNodeCount;
+ Dimension_Index dim = forceVector->dim;
+ FeMesh* mesh = forceVector->feVariable->feMesh;
+ Grid* grid;
+ Node_ElementLocalIndex eNode_I;
+ ElementType* elementType;
+ Dof_Index nodeDofCount;
+ double stress, area;
+ IJK ijk;
+
+ Node_DomainIndex* elementNodes = NULL;
+
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ nodeDofCount = dim;
+
+ grid = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+
+ area=StressBC_compute_face_area(self->_wall,mesh,lElement_I,
+ &elementNodeCount,dim,&elementNodes);
+ /* Apply the stress */
+ for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
+ /* Make sure that we are on the boundary */
+ int condition, entry_I;
+ ConditionFunction* cf;
+ RegularMeshUtils_Node_1DTo3D
+ ( mesh, Mesh_DomainToGlobal(mesh, MT_VERTEX,
+ elementNodes[eNode_I]), ijk );
+ switch(self->_wall)
+ {
+ case Wall_Left:
+ condition=(ijk[0] == 0);
+ break;
+ case Wall_Right:
+ condition=(ijk[0] == ( grid->sizes[0] - 1 ));
+ break;
+ case Wall_Bottom:
+ condition=(ijk[1] == 0);
+ break;
+ case Wall_Top:
+ condition=(ijk[1] == ( grid->sizes[1] - 1 ));
+ break;
+ case Wall_Front:
+ condition=(ijk[2] == 0);
+ break;
+ case Wall_Back:
+ condition=(ijk[2] == ( grid->sizes[2] - 1 ));
+ break;
+ }
+
+ if(condition)
+ {
+ unsigned overcount=StressBC_get_overcount(dim,ijk,grid->sizes);
+ for(entry_I=0; entry_I<self->numEntries; ++entry_I)
+ {
+ switch(self->_entryTbl[entry_I].type)
+ {
+ case StressBC_Double:
+ stress=self->_entryTbl[entry_I].DoubleValue;
+ break;
+ case StressBC_ConditionFunction:
+ cf=self->conFunc_Register->
+ _cf[self->_entryTbl[entry_I].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,&stress);
+ break;
+ }
+ /* We have to divide by an overcount_factor, because
+ otherwise different elements will count the same
+ node more than once. In 2D, nodes are counted
+ twice, in 3D, nodes are counted four times. */
+ elForceVec[ eNode_I * nodeDofCount
+ + self->_entryTbl[entry_I].axis ] += stress*area/overcount;
+ }
+ }
+ }
+}
+
+
+double StressBC_compute_face_area(Wall wall, FeMesh *mesh, Index lElement_I,
+ Element_NodeIndex *elementNodeCount,
+ Dimension_Index dim,
+ Node_DomainIndex **elementNodes)
+{
+ /* Compute the area of the face. */
+ if(dim==2)
+ {
+ double *coord1, *coord2;
+ int lower,upper,direction;
+ switch(wall)
+ {
+ case Wall_Left:
+ lower=0;
+ upper=3;
+ direction=1;
+ break;
+ case Wall_Right:
+ lower=1;
+ upper=2;
+ direction=1;
+ break;
+ case Wall_Bottom:
+ lower=0;
+ upper=1;
+ direction=0;
+ break;
+ case Wall_Top:
+ lower=3;
+ upper=2;
+ direction=0;
+ break;
+ }
+
+ Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
+ MT_VERTEX,elementNodeCount, elementNodes);
+
+ coord1=Mesh_GetVertex(mesh,(*elementNodes)[lower]);
+ coord2=Mesh_GetVertex(mesh,(*elementNodes)[upper]);
+ return coord2[direction]-coord1[direction];
+ }
+ else
+ {
+ double *coord1, *coord2, *coord3, *coord4;
+ int first, second, third, fourth;
+
+ /* StGermain uses the ordering
+ 0: x,y,z
+ 1: x+,y,z
+ 2: x,y+,z
+ 3: x+,y+,z
+ 4: x,y,z+
+ 5: x+,y,z+
+ 6: x,y+,z+
+ 7: x+,y+,z+
+ */
+
+ /* Get the indices for which wall we want to get the area of. */
+ switch(wall)
+ {
+ case Wall_Left:
+ first=0;
+ second=2;
+ third=6;
+ fourth=4;
+ break;
+ case Wall_Right:
+ first=1;
+ second=3;
+ third=7;
+ fourth=5;
+ break;
+ case Wall_Bottom:
+ first=0;
+ second=1;
+ third=5;
+ fourth=4;
+ break;
+ case Wall_Top:
+ first=2;
+ second=3;
+ third=7;
+ fourth=6;
+ break;
+ case Wall_Front:
+ first=0;
+ second=1;
+ third=4;
+ fourth=3;
+ break;
+ case Wall_Back:
+ first=4;
+ second=5;
+ third=7;
+ fourth=6;
+ break;
+ }
+
+ Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
+ MT_VERTEX,elementNodeCount, elementNodes);
+
+ coord1=Mesh_GetVertex(mesh,(*elementNodes)[first]);
+ coord2=Mesh_GetVertex(mesh,(*elementNodes)[second]);
+ coord3=Mesh_GetVertex(mesh,(*elementNodes)[third]);
+ coord4=Mesh_GetVertex(mesh,(*elementNodes)[fourth]);
+
+ return StGermain_ConvexQuadrilateralArea(coord1,coord2,coord3,coord4,
+ dim);
+ }
+}
+
+/* We have to divide by an overcount_factor, because otherwise
+ different elements will count the same node more than once. In 2D,
+ in the interior, there are 4 elements that touch a node. On an
+ edge, 2, and on the corner 1. Similary for 3D. */
+unsigned StressBC_get_overcount(Dimension_Index dim, IJK ijk, unsigned sizes[])
+{
+ unsigned overcount, i;
+ if(dim==2)
+ overcount=4;
+ else
+ overcount=8;
+
+ for(i=0;i<dim;++i)
+ if(ijk[i]==0 || ijk[i]==sizes[dim]-1)
+ overcount/=2;
+ return overcount;
+}
Added: long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.h 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.h 2007-07-18 18:37:29 UTC (rev 7696)
@@ -0,0 +1,117 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+** Melbourne, 3053, Australia.
+** Copyright (c) 2005-2006, Monash Cluster Computing, Building 28, Monash University Clayton Campus,
+** Victoria, 3800, Australia
+**
+** Primary Contributing Organisations:
+** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+** Australian Computational Earth Systems Simulator - http://www.access.edu.au
+** Monash Cluster Computing - http://www.mcc.monash.edu.au
+**
+** Contributors:
+** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+** Patrick D. Sunter, Software Engineer, VPAC. (patrick at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** David Stegman, Postdoctoral Fellow, Monash University. (david.stegman at sci.monash.edu.au)
+** Wendy Sharples, PhD Student, Monash University (wendy.sharples at sci.monash.edu.au)
+**
+** Modified 2006 by Walter Landry (California Institute of Technology)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+*/
+
+
+#ifndef __Gale_Utils_StressBC_h__
+#define __Gale_Utils_StressBC_h__
+
+ /** Textual name of this class */
+ extern const Type StressBC_Type;
+
+ /** StressBC class contents */
+ #define __StressBC \
+ /* General info */ \
+ __ForceTerm \
+ \
+ /* StressBC info */ \
+ Wall _wall; \
+ StressBC_Entry _entryTbl[3]; \
+ int numEntries; \
+ ConditionFunction_Register* conFunc_Register; \
+ FiniteElementContext* context; \
+
+ struct StressBC { __StressBC };
+
+ StressBC* StressBC_New(
+ Name name,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context);
+
+ StressBC* _StressBC_New(
+ SizeT sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ ForceTerm_AssembleElementFunction* _assembleElement,
+ ConditionFunction_Register* conFunc_Register,
+ Name name,
+ FiniteElementContext* context);
+
+ void StressBC_InitAll(
+ void* forceTerm,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context);
+
+ void _StressBC_Delete( void* forceTerm );
+ void _StressBC_Print( void* forceTerm, Stream* stream );
+
+ void* _StressBC_DefaultNew( Name name ) ;
+ void _StressBC_Construct( void* forceTerm, Stg_ComponentFactory* cf, void* data ) ;
+ void _StressBC_Build( void* forceTerm, void* data ) ;
+ void _StressBC_Initialise( void* forceTerm, void* data ) ;
+ void _StressBC_Execute( void* forceTerm, void* data ) ;
+ void _StressBC_Destroy( void* forceTerm, void* data ) ;
+
+ void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) ;
+ void _StressBC_GetValues(Stg_ComponentFactory* cf, void *stressBC, char *direction);
+ double StressBC_compute_face_area(Wall wall, FeMesh *mesh,
+ Index lElement_I,
+ Element_NodeIndex *elementNodeCount,
+ Dimension_Index dim,
+ Node_DomainIndex **elementNodes);
+ unsigned StressBC_get_overcount(Dimension_Index dim, IJK ijk,
+ unsigned sizes[]);
+
+#endif
Added: long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.meta
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.meta 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.meta 2007-07-18 18:37:29 UTC (rev 7696)
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">StressBC</param>
+<param name="Organisation">VPAC and MCC</param>
+<param name="Project">Gale</param>
+<param name="Location">./Gale/Utils/src/</param>
+<param name="Project Web">https://csd.vpac.org/twiki/bin/view/PICellerator/WebHome</param>
+<param name="Copyright">Copyright (C) 2005 VPAC and Monash Cluster Computing.</param>
+<param name="License">https://csd.vpac.org/twiki/bin/view/Stgermain/SoftwareLicense http://www.opensource.org/licenses/bsd-license.php</param>
+<param name="Parent">ForceTerm</param>
+<param name="Description">This adds a stress boundary condition.</param>
+
+<!--Now the interesting stuff-->
+
+
+<list name="Params">
+ <struct>
+ </struct>
+
+</list>
+
+<!-- Add an exmaple XML if possible -->
+<param name="Example">
+ <struct name="stressBC">
+ <param name="Type">StressBC</param>
+ <param name="ForceVector">mom_force</param>
+ <param name="Swarm">materialPoints</param>
+ <param name="wall">bottom</param>
+ <param name="x_type">double</param>
+ <param name="x_value">1.0</param>
+ </struct>
+</param>
+
+</StGermainData>
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/Utils.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/Utils.h 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/Utils.h 2007-07-18 18:37:29 UTC (rev 7696)
@@ -51,5 +51,7 @@
#include "Init.h"
#include "Finalise.h"
#include "FrictionVC.h"
+ #include "KineticFriction.h"
+ #include "StressBC.h"
#endif
Added: long/3D/Gale/trunk/src/Gale/Utils/src/Wall.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/Wall.c 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/Wall.c 2007-07-18 18:37:29 UTC (rev 7696)
@@ -0,0 +1,12 @@
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StgFEM/StgFEM.h>
+#include "types.h"
+
+const char* WallEnumToStr[Wall_Size] = {
+ "back",
+ "left",
+ "bottom",
+ "right",
+ "top",
+ "front" };
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/types.h 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/types.h 2007-07-18 18:37:29 UTC (rev 7696)
@@ -46,22 +46,38 @@
#ifndef __Gale_Utils_types_h__
#define __Gale_Utils_types_h__
- typedef struct GaleContext GaleContext;
- typedef struct _FrictionVC FrictionVC;
+typedef struct GaleContext GaleContext;
+typedef struct _FrictionVC FrictionVC;
+typedef struct KineticFriction KineticFriction;
+typedef struct StressBC StressBC;
- /* FrictionVC types */
- typedef enum
- {
- FrictionVC_Wall_Back,
- FrictionVC_Wall_Left,
- FrictionVC_Wall_Bottom,
- FrictionVC_Wall_Right,
- FrictionVC_Wall_Top,
- FrictionVC_Wall_Front,
- FrictionVC_Wall_Size
- } FrictionVC_Wall;
-
- typedef struct _FrictionVC_Entry FrictionVC_Entry;
- typedef Index FrictionVC_Entry_Index;
+typedef enum
+ {
+ StressBC_Double,
+ StressBC_ConditionFunction
+ } StressBC_Types;
+typedef struct
+{
+ StressBC_Types type;
+ double DoubleValue;
+ Index CFIndex;
+ Axis axis;
+} StressBC_Entry;
+
+/* Wall types */
+typedef enum
+ {
+ Wall_Back,
+ Wall_Left,
+ Wall_Bottom,
+ Wall_Right,
+ Wall_Top,
+ Wall_Front,
+ Wall_Size
+ } Wall;
+
+typedef struct _FrictionVC_Entry FrictionVC_Entry;
+typedef Index FrictionVC_Entry_Index;
+
#endif
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/Init.c
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/Init.c 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/Init.c 2007-07-18 18:37:29 UTC (rev 7696)
@@ -64,9 +64,6 @@
RegisterParent( BuoyancyForceTerm_Type, ForceTerm_Type );
RegisterParent( BuoyancyForceTermThermoChem_Type, ForceTerm_Type );
- Stg_ComponentRegister_Add( componentsRegister, StressBC_Type, "0", _StressBC_DefaultNew );
-
- RegisterParent( StressBC_Type, ForceTerm_Type );
RegisterParent( MaterialSwarmVariable_Type, SwarmVariable_Type );
RegisterParent( PCDVC_Type, DVCWeights_Type );
return True;
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/SConscript
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/SConscript 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/SConscript 2007-07-18 18:37:29 UTC (rev 7696)
@@ -25,7 +25,6 @@
Init.h
MaterialSwarmVariable.h
PCDVC.h
-StressBC.h
types.h
Utils.h""")
@@ -35,14 +34,12 @@
Finalise.c
Init.c
PCDVC.c
-StressBC.c
MaterialSwarmVariable.c""")
meta_files=Split("""BuoyancyForceTerm.meta
BuoyancyForceTermThermoChem.meta
PCDVC.meta
-StressBC.meta
MaterialSwarmVariable.meta""")
c_files+=[local_env.Meta(meta_files)]
Deleted: long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c 2007-07-18 18:37:29 UTC (rev 7696)
@@ -1,540 +0,0 @@
-/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-**
-** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
-** Melbourne, 3053, Australia.
-** Copyright (c) 2005-2006, Monash Cluster Computing, Building 28, Monash University Clayton Campus,
-** Victoria, 3800, Australia
-**
-** Primary Contributing Organisations:
-** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
-** Australian Computational Earth Systems Simulator - http://www.access.edu.au
-** Monash Cluster Computing - http://www.mcc.monash.edu.au
-**
-** Contributors:
-** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
-** Patrick D. Sunter, Software Engineer, VPAC. (patrick at vpac.org)
-** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
-** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
-** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
-** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
-** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
-** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
-** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
-** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
-** David Stegman, Postdoctoral Fellow, Monash University. (david.stegman at sci.monash.edu.au)
-** Wendy Sharples, PhD Student, Monash University (wendy.sharples at sci.monash.edu.au)
-**
-** Modified 2006 by Walter Landry (California Institute of Technology)
-**
-** This library is free software; you can redistribute it and/or
-** modify it under the terms of the GNU Lesser General Public
-** License as published by the Free Software Foundation; either
-** version 2.1 of the License, or (at your option) any later version.
-**
-** This library is distributed in the hope that it will be useful,
-** but WITHOUT ANY WARRANTY; without even the implied warranty of
-** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-** Lesser General Public License for more details.
-**
-** You should have received a copy of the GNU Lesser General Public
-** License along with this library; if not, write to the Free Software
-** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-**
-** $Id$
-**
-**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
-
-#include <mpi.h>
-#include <StGermain/StGermain.h>
-#include <StgFEM/StgFEM.h>
-#include <PICellerator/Voronoi/Voronoi.h>
-#include <PICellerator/PopulationControl/PopulationControl.h>
-#include <PICellerator/Weights/Weights.h>
-#include <PICellerator/MaterialPoints/MaterialPoints.h>
-
-#include "PICellerator/Utils/types.h"
-#include "PICellerator/Utils/MaterialSwarmVariable.h"
-
-#include "types.h"
-#include "StressBC.h"
-
-#include <assert.h>
-#include <string.h>
-
-/* Textual name of this class */
-const Type StressBC_Type = "StressBC";
-
-const char* StressBC_WallEnumToStr[StressBC_Wall_Size] = {
- "back",
- "left",
- "bottom",
- "right",
- "top",
- "front" };
-
-StressBC* StressBC_New(
- Name name,
- ForceVector* forceVector,
- Swarm* integrationSwarm,
- ConditionFunction_Register* conFunc_Register,
- FiniteElementContext* context)
-{
- StressBC* self = (StressBC*) _StressBC_DefaultNew( name );
-
- StressBC_InitAll(
- self,
- forceVector,
- integrationSwarm,
- conFunc_Register,
- context);
-
- return self;
-}
-
-/* Creation implementation / Virtual constructor */
-StressBC* _StressBC_New(
- SizeT sizeOfSelf,
- Type type,
- Stg_Class_DeleteFunction* _delete,
- Stg_Class_PrintFunction* _print,
- Stg_Class_CopyFunction* _copy,
- Stg_Component_DefaultConstructorFunction* _defaultConstructor,
- Stg_Component_ConstructFunction* _construct,
- Stg_Component_BuildFunction* _build,
- Stg_Component_InitialiseFunction* _initialise,
- Stg_Component_ExecuteFunction* _execute,
- Stg_Component_DestroyFunction* _destroy,
- ForceTerm_AssembleElementFunction* _assembleElement,
- ConditionFunction_Register* conFunc_Register,
- Name name,
- FiniteElementContext* context)
-{
- StressBC* self;
-
- /* Allocate memory */
- assert( sizeOfSelf >= sizeof(StressBC) );
- self = (StressBC*) _ForceTerm_New(
- sizeOfSelf,
- type,
- _delete,
- _print,
- _copy,
- _defaultConstructor,
- _construct,
- _build,
- _initialise,
- _execute,
- _destroy,
- _assembleElement,
- name );
-
- self->conFunc_Register=conFunc_Register;
- self->context=context;
-
- return self;
-}
-
-void _StressBC_Init(
- StressBC* self,
- ConditionFunction_Register* conFunc_Register,
- FiniteElementContext* context)
-{
- self->numEntries = 0;
- self->conFunc_Register=conFunc_Register;
- self->context=context;
-}
-
-void StressBC_InitAll(
- void* forceTerm,
- ForceVector* forceVector,
- Swarm* integrationSwarm,
- ConditionFunction_Register* conFunc_Register,
- FiniteElementContext* context)
-{
- StressBC* self = (StressBC*) forceTerm;
-
- ForceTerm_InitAll( self, forceVector, integrationSwarm, NULL );
- _StressBC_Init( self, conFunc_Register, context );
-}
-
-void _StressBC_Delete( void* forceTerm ) {
- StressBC* self = (StressBC*)forceTerm;
-
- _ForceTerm_Delete( self );
-}
-
-void _StressBC_Print( void* forceTerm, Stream* stream ) {
- StressBC* self = (StressBC*)forceTerm;
-
- _ForceTerm_Print( self, stream );
-}
-
-void* _StressBC_DefaultNew( Name name ) {
- return (void*)_StressBC_New(
- sizeof(StressBC),
- StressBC_Type,
- _StressBC_Delete,
- _StressBC_Print,
- NULL,
- _StressBC_DefaultNew,
- _StressBC_Construct,
- _StressBC_Build,
- _StressBC_Initialise,
- _StressBC_Execute,
- _StressBC_Destroy,
- _StressBC_AssembleElement,
- NULL,
- name,
- NULL);
-}
-
-void _StressBC_Construct( void* forceTerm, Stg_ComponentFactory* cf, void* data ) {
- StressBC* self = (StressBC*)forceTerm;
- FiniteElementContext* context;
-
- /* Construct Parent */
- _ForceTerm_Construct( self, cf, data );
-
- context = (FiniteElementContext*)Stg_ComponentFactory_ConstructByName
- ( cf, "context", FiniteElementContext, True, data ) ;
-
- _StressBC_Init( self, context->condFunc_Register, context );
- {
- char* wallStr;
-
- /* Obtain which wall */
- wallStr = Stg_ComponentFactory_GetString( cf, self->name,
- "wall", "");
- if (!strcasecmp(wallStr, "back"))
- self->_wall = StressBC_Wall_Back;
- else if (!strcasecmp(wallStr, "left"))
- self->_wall = StressBC_Wall_Left;
- else if (!strcasecmp(wallStr, "bottom"))
- self->_wall = StressBC_Wall_Bottom;
- else if (!strcasecmp(wallStr, "right"))
- self->_wall = StressBC_Wall_Right;
- else if (!strcasecmp(wallStr, "top"))
- self->_wall = StressBC_Wall_Top;
- else if (!strcasecmp(wallStr, "front"))
- self->_wall = StressBC_Wall_Front;
- else {
- assert( 0 );
- self->_wall = StressBC_Wall_Size; /* invalid entry */
- }
- }
- _StressBC_GetValues(cf,self,"x");
- _StressBC_GetValues(cf,self,"y");
- _StressBC_GetValues(cf,self,"z");
-
-}
-
-/* Gets the actual values used by the StressBC (e.g. a float or a function). */
-void _StressBC_GetValues(Stg_ComponentFactory* cf, void *stressBC,
- char *direction)
-{
- StressBC* self = (StressBC*)stressBC;
- char temp_str[20];
- char *type;
-
- switch(*direction)
- {
- case 'x':
- self->_entryTbl[self->numEntries].axis=I_AXIS;
- break;
- case 'y':
- self->_entryTbl[self->numEntries].axis=J_AXIS;
- break;
- case 'z':
- self->_entryTbl[self->numEntries].axis=K_AXIS;
- break;
- }
-
- strcat(strcpy(temp_str,direction),"_type");
- type=Stg_ComponentFactory_GetString( cf, self->name, temp_str, "");
- strcat(strcpy(temp_str,direction),"_value");
-
- if(!strcasecmp(type,"double") || !strcasecmp(type,"float"))
- {
- self->_entryTbl[self->numEntries].type = StressBC_Double;
- self->_entryTbl[self->numEntries].DoubleValue =
- Stg_ComponentFactory_GetDouble( cf, self->name, temp_str, 0.0);
- (self->numEntries)++;
- }
- else if(!strcasecmp(type,"func"))
- {
- char *funcName =
- Stg_ComponentFactory_GetString( cf, self->name, temp_str, "");
- Index cfIndex;
-
- cfIndex = ConditionFunction_Register_GetIndex
- ( self->conFunc_Register, funcName);
- self->_entryTbl[self->numEntries].type =
- StressBC_ConditionFunction;
- if ( cfIndex == (unsigned)-1 ) {
- Stream* errorStr = Journal_Register( Error_Type, self->type );
-
- Journal_Printf( errorStr, "Error- in %s: While parsing "
- "definition of StressBC (applies to wall \"%s\"), the cond. func. applied to "
- "direction \"%s\" - \"%s\" - wasn't found in the c.f. register.\n",
- __func__, StressBC_WallEnumToStr[self->_wall],
- direction, funcName );
- Journal_Printf( errorStr, "(Available functions in the C.F. register are: ");
- ConditionFunction_Register_PrintNameOfEachFunc
- ( self->conFunc_Register, errorStr );
- Journal_Printf( errorStr, ")\n");
- assert(0);
- }
- self->_entryTbl[self->numEntries].CFIndex = cfIndex;
- (self->numEntries)++;
- }
- else if(strlen(type)!=0)
- {
- Stream* errorStr = Journal_Register( Error_Type, self->type );
- Journal_Printf( errorStr, "Error- in %s: While parsing "
- "definition of StressBC (applies to wall \"%s\"), the type of condition \"%s\"\nis not supported. Supported types are \"double\" and \"function\".",
- __func__, StressBC_WallEnumToStr[self->_wall],
- type );
- assert(0);
- }
-}
-
-void _StressBC_Build( void* forceTerm, void* data ) {
- StressBC* self = (StressBC*)forceTerm;
- Name name;
-
- _ForceTerm_Build( self, data );
-}
-
-void _StressBC_Initialise( void* forceTerm, void* data ) {
- StressBC* self = (StressBC*)forceTerm;
- Index i;
-
- _ForceTerm_Initialise( self, data );
-
-}
-
-void _StressBC_Execute( void* forceTerm, void* data ) {
- _ForceTerm_Execute( forceTerm, data );
-}
-
-void _StressBC_Destroy( void* forceTerm, void* data ) {
- _ForceTerm_Destroy( forceTerm, data );
-}
-
-double cross(double *coord1, double *coord2, int normal)
-{
- switch(normal)
- {
- case 0:
- return coord1[1]*coord2[2]-coord1[2]*coord2[1];
- break;
- case 1:
- return coord1[0]*coord2[2]-coord1[2]*coord2[0];
- break;
- case 2:
- return coord1[0]*coord2[1]-coord1[1]*coord2[0];
- break;
- default:
- printf("Bad call of cross() in StressBC\n");
- abort();
- }
- return 0;
-}
-
-void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
- StressBC* self = (StressBC*) forceTerm;
- Element_NodeIndex elementNodeCount;
- Dimension_Index dim = forceVector->dim;
- FeMesh* mesh = forceVector->feVariable->feMesh;
- Grid* grid;
- Node_ElementLocalIndex eNode_I;
- ElementType* elementType;
- Dof_Index nodeDofCount;
- double stress, area;
- IJK ijk;
- int overcount_factor;
-
- Node_DomainIndex* elementNodes = NULL;
-
- elementType = FeMesh_GetElementType( mesh, lElement_I );
- nodeDofCount = dim;
-
- grid = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
- ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
-
- /* Compute the area of the face. */
- if(dim==2)
- {
- double *coord1, *coord2;
- int lower,upper,direction;
- switch(self->_wall)
- {
- case StressBC_Wall_Left:
- lower=0;
- upper=3;
- direction=1;
- break;
- case StressBC_Wall_Right:
- lower=1;
- upper=2;
- direction=1;
- break;
- case StressBC_Wall_Bottom:
- lower=0;
- upper=1;
- direction=0;
- break;
- case StressBC_Wall_Top:
- lower=3;
- upper=2;
- direction=0;
- break;
- }
-
- Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
- MT_VERTEX,&elementNodeCount, &elementNodes);
-
- coord1=Mesh_GetVertex(mesh,elementNodes[lower]);
- coord2=Mesh_GetVertex(mesh,elementNodes[upper]);
- area=coord2[direction]-coord1[direction];
- overcount_factor=2;
- }
- else
- {
- double *coord1, *coord2, *coord3, *coord4;
- int first, second, third, fourth, normal;
-
-
- Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
- MT_VERTEX,&elementNodeCount, &elementNodes);
-
- /* StGermain uses the ordering
- 0: x,y,z
- 1: x+,y,z
- 2: x,y+,z
- 3: x+,y+,z
- 4: x,y,z+
- 5: x+,y,z+
- 6: x,y+,z+
- 7: x+,y+,z+
- */
-
- switch(self->_wall)
- {
- case StressBC_Wall_Left:
- first=0;
- second=2;
- third=6;
- fourth=4;
- normal=0;
- break;
- case StressBC_Wall_Right:
- first=1;
- second=3;
- third=7;
- fourth=5;
- normal=0;
- break;
- case StressBC_Wall_Bottom:
- first=0;
- second=1;
- third=5;
- fourth=4;
- normal=1;
- break;
- case StressBC_Wall_Top:
- first=2;
- second=3;
- third=7;
- fourth=6;
- normal=1;
- break;
- case StressBC_Wall_Front:
- first=0;
- second=1;
- third=4;
- fourth=3;
- normal=2;
- break;
- case StressBC_Wall_Back:
- first=4;
- second=5;
- third=7;
- fourth=6;
- normal=2;
- break;
- }
-
- Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
- MT_VERTEX,&elementNodeCount, &elementNodes);
-
- coord1=Mesh_GetVertex(mesh,elementNodes[first]);
- coord2=Mesh_GetVertex(mesh,elementNodes[second]);
- coord3=Mesh_GetVertex(mesh,elementNodes[third]);
- coord4=Mesh_GetVertex(mesh,elementNodes[fourth]);
- area=fabs(cross(coord1,coord2,normal) + cross(coord2,coord3,normal)
- + cross(coord3,coord4,normal)
- + cross(coord4,coord1,normal))/2;
- overcount_factor=4;
- }
-
- /* Apply the stress */
- for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
- /* Make sure that we are on the boundary */
- int condition, entry_I;
- ConditionFunction* cf;
- RegularMeshUtils_Node_1DTo3D( mesh, Mesh_DomainToGlobal(mesh, MT_VERTEX, elementNodes[eNode_I]), ijk );
- switch(self->_wall)
- {
- case StressBC_Wall_Left:
- condition=(ijk[0] == 0);
- break;
- case StressBC_Wall_Right:
- condition=(ijk[0] == ( grid->sizes[0] - 1 ));
- break;
- case StressBC_Wall_Bottom:
- condition=(ijk[1] == 0);
- break;
- case StressBC_Wall_Top:
- condition=(ijk[1] == ( grid->sizes[1] - 1 ));
- break;
- case StressBC_Wall_Front:
- condition=(ijk[2] == 0);
- break;
- case StressBC_Wall_Back:
- condition=(ijk[2] == ( grid->sizes[2] - 1 ));
- break;
- }
-
- if(condition)
- {
- for(entry_I=0; entry_I<self->numEntries; ++entry_I)
- {
- switch(self->_entryTbl[entry_I].type)
- {
- case StressBC_Double:
- stress=self->_entryTbl[entry_I].DoubleValue;
- break;
- case StressBC_ConditionFunction:
- cf=self->conFunc_Register->
- _cf[self->_entryTbl[entry_I].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,&stress);
- break;
- }
- /* We have to divide by an overcount_factor, because
- otherwise different elements will count the same
- node more than once. In 2D, nodes are counted
- twice, in 3D, nodes are counted four times. */
- elForceVec[ eNode_I * nodeDofCount
- + self->_entryTbl[entry_I].axis ] += stress*area/overcount_factor;
- }
- }
- }
-}
-
Deleted: long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.h
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.h 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.h 2007-07-18 18:37:29 UTC (rev 7696)
@@ -1,110 +0,0 @@
-/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-**
-** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
-** Melbourne, 3053, Australia.
-** Copyright (c) 2005-2006, Monash Cluster Computing, Building 28, Monash University Clayton Campus,
-** Victoria, 3800, Australia
-**
-** Primary Contributing Organisations:
-** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
-** Australian Computational Earth Systems Simulator - http://www.access.edu.au
-** Monash Cluster Computing - http://www.mcc.monash.edu.au
-**
-** Contributors:
-** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
-** Patrick D. Sunter, Software Engineer, VPAC. (patrick at vpac.org)
-** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
-** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
-** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
-** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
-** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
-** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
-** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
-** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
-** David Stegman, Postdoctoral Fellow, Monash University. (david.stegman at sci.monash.edu.au)
-** Wendy Sharples, PhD Student, Monash University (wendy.sharples at sci.monash.edu.au)
-**
-** Modified 2006 by Walter Landry (California Institute of Technology)
-**
-** This library is free software; you can redistribute it and/or
-** modify it under the terms of the GNU Lesser General Public
-** License as published by the Free Software Foundation; either
-** version 2.1 of the License, or (at your option) any later version.
-**
-** This library is distributed in the hope that it will be useful,
-** but WITHOUT ANY WARRANTY; without even the implied warranty of
-** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-** Lesser General Public License for more details.
-**
-** You should have received a copy of the GNU Lesser General Public
-** License along with this library; if not, write to the Free Software
-** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-**
-*/
-
-
-#ifndef __Gale_Utils_StressBC_h__
-#define __Gale_Utils_StressBC_h__
-
- /** Textual name of this class */
- extern const Type StressBC_Type;
-
- /** StressBC class contents */
- #define __StressBC \
- /* General info */ \
- __ForceTerm \
- \
- /* StressBC info */ \
- StressBC_Wall _wall; \
- StressBC_Entry _entryTbl[3]; \
- int numEntries; \
- ConditionFunction_Register* conFunc_Register; \
- FiniteElementContext* context; \
-
- struct StressBC { __StressBC };
-
- StressBC* StressBC_New(
- Name name,
- ForceVector* forceVector,
- Swarm* integrationSwarm,
- ConditionFunction_Register* conFunc_Register,
- FiniteElementContext* context);
-
- StressBC* _StressBC_New(
- SizeT sizeOfSelf,
- Type type,
- Stg_Class_DeleteFunction* _delete,
- Stg_Class_PrintFunction* _print,
- Stg_Class_CopyFunction* _copy,
- Stg_Component_DefaultConstructorFunction* _defaultConstructor,
- Stg_Component_ConstructFunction* _construct,
- Stg_Component_BuildFunction* _build,
- Stg_Component_InitialiseFunction* _initialise,
- Stg_Component_ExecuteFunction* _execute,
- Stg_Component_DestroyFunction* _destroy,
- ForceTerm_AssembleElementFunction* _assembleElement,
- ConditionFunction_Register* conFunc_Register,
- Name name,
- FiniteElementContext* context);
-
- void StressBC_InitAll(
- void* forceTerm,
- ForceVector* forceVector,
- Swarm* integrationSwarm,
- ConditionFunction_Register* conFunc_Register,
- FiniteElementContext* context);
-
- void _StressBC_Delete( void* forceTerm );
- void _StressBC_Print( void* forceTerm, Stream* stream );
-
- void* _StressBC_DefaultNew( Name name ) ;
- void _StressBC_Construct( void* forceTerm, Stg_ComponentFactory* cf, void* data ) ;
- void _StressBC_Build( void* forceTerm, void* data ) ;
- void _StressBC_Initialise( void* forceTerm, void* data ) ;
- void _StressBC_Execute( void* forceTerm, void* data ) ;
- void _StressBC_Destroy( void* forceTerm, void* data ) ;
-
- void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) ;
- void _StressBC_GetValues(Stg_ComponentFactory* cf, void *stressBC, char *direction);
-
-#endif
Deleted: long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.meta
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.meta 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.meta 2007-07-18 18:37:29 UTC (rev 7696)
@@ -1,36 +0,0 @@
-<?xml version="1.0"?>
-<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
-<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
-
-<param name="Name">StressBC</param>
-<param name="Organisation">VPAC and MCC</param>
-<param name="Project">Gale</param>
-<param name="Location">./Gale/Utils/src/</param>
-<param name="Project Web">https://csd.vpac.org/twiki/bin/view/PICellerator/WebHome</param>
-<param name="Copyright">Copyright (C) 2005 VPAC and Monash Cluster Computing.</param>
-<param name="License">https://csd.vpac.org/twiki/bin/view/Stgermain/SoftwareLicense http://www.opensource.org/licenses/bsd-license.php</param>
-<param name="Parent">ForceTerm</param>
-<param name="Description">This adds a stress boundary condition.</param>
-
-<!--Now the interesting stuff-->
-
-
-<list name="Params">
- <struct>
- </struct>
-
-</list>
-
-<!-- Add an exmaple XML if possible -->
-<param name="Example">
- <struct name="stressBC">
- <param name="Type">StressBC</param>
- <param name="ForceVector">mom_force</param>
- <param name="Swarm">materialPoints</param>
- <param name="wall">bottom</param>
- <param name="x_type">double</param>
- <param name="x_value">1.0</param>
- </struct>
-</param>
-
-</StGermainData>
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/Utils.h
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/Utils.h 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/Utils.h 2007-07-18 18:37:29 UTC (rev 7696)
@@ -46,7 +46,6 @@
#include "types.h"
#include "BuoyancyForceTerm.h"
- #include "StressBC.h"
#include "BuoyancyForceTermThermoChem.h"
#include "MaterialSwarmVariable.h"
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/types.h 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/types.h 2007-07-18 18:37:29 UTC (rev 7696)
@@ -55,35 +55,8 @@
#define __PICellerator_Utils_types_h__
typedef struct BuoyancyForceTerm BuoyancyForceTerm;
- typedef struct StressBC StressBC;
typedef struct MaterialSwarmVariable MaterialSwarmVariable;
typedef struct BuoyancyForceTermThermoChem BuoyancyForceTermThermoChem;
typedef struct PCDVC PCDVC;
-/* StressBC types */
-typedef enum
- {
- StressBC_Wall_Back,
- StressBC_Wall_Left,
- StressBC_Wall_Bottom,
- StressBC_Wall_Right,
- StressBC_Wall_Top,
- StressBC_Wall_Front,
- StressBC_Wall_Size
- } StressBC_Wall;
-
-typedef enum
- {
- StressBC_Double,
- StressBC_ConditionFunction
- } StressBC_Types;
-
-typedef struct
-{
- StressBC_Types type;
- double DoubleValue;
- Index CFIndex;
- Axis axis;
-} StressBC_Entry;
-
#endif
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c 2007-07-18 18:37:22 UTC (rev 7695)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c 2007-07-18 18:37:29 UTC (rev 7696)
@@ -669,8 +669,11 @@
/* Add this to create a new set of bcs dependent on the stress used in friction bc's. */
stress = (FeVariable*)FieldVariable_Register_GetByName( ((FiniteElementContext *)_context)->fieldVariable_Register, "StressField" );
- ParticleFeVariable_Update( stress );
- _SystemLinearEquations_Build(sle,_context);
+ if(stress!=NULL)
+ {
+ ParticleFeVariable_Update( stress );
+ _SystemLinearEquations_Build(sle,_context);
+ }
self->linearExecute( self, _context );
self->hasExecuted = True;
@@ -683,8 +686,11 @@
for ( self->nonLinearIteration_I = 1 ; self->nonLinearIteration_I < maxIterations ; self->nonLinearIteration_I++ ) {
/* Add this to create a new set of bcs dependent on the stress used in friction bc's. */
- ParticleFeVariable_Update( stress );
- _SystemLinearEquations_Build(sle,_context);
+ if(stress!=NULL)
+ {
+ ParticleFeVariable_Update( stress );
+ _SystemLinearEquations_Build(sle,_context);
+ }
Vector_CopyEntries( currentVector, previousVector );
More information about the cig-commits
mailing list