sandbox/bugs/tensor_boundary.c
Let’s describe a 2D tensor as: (a) a tensor, (b) 2 vectors or (c) 4 scalars
tensor S[];
vector TVx[], TVy[];
scalar TSxx[], TSxy[], TSyx[], TSyy[];
Let’s apply only a shear in x direction (\tau_{xy} = 1). The rest of components of the tensor are zero.
int main()
{
init_grid (4);
vector v = S.x;
v.n[top] = dirichlet(1); // Sxy
v.t[top] = dirichlet(0.); // Sxx
v = S.y;
v.n[top] = dirichlet(0); // Syy
v.t[top] = dirichlet(0); // Syx
TVx.n[top] = dirichlet(1);
TVx.t[top] = dirichlet(0);
TVy.n[top] = dirichlet(0);
TVy.t[top] = dirichlet(0);
TSxx[top] = dirichlet(0);
TSxy[top] = dirichlet(1);
TSyx[top] = dirichlet(0);
TSyy[top] = dirichlet(0);
foreach() {
foreach_dimension () {
S.x.x[] = 0;
S.x.y[] = 0;
TVx.x[] = 0.;
TVy.x[] = 0.;
}
TSxx[] = 0;
TSxy[] = 0;
TSyx[] = 0;
TSyy[] = 0;
}
boundary ((scalar *){S.x, S.y, TVx, TVy, TSxx, TSxy, TSyx, TSyy});
double shear, normal;
Lets compute the x component of the corresponding stress (acceleration): \displaystyle \nabla \cdot \tau |_x =(\tau_xx)_x + (\tau_{xy})_y The shear component should be the same despite the tensor was constructed in (a), (b) or (c) manner.
foreach_face(x) {
shear = (S.x.y[0,1]+S.x.y[-1,1]-S.x.y[0,-1]-S.x.y[-1,-1])/(4.*Delta); //(\tau_{xy})_y
normal = (S.x.x[]-S.x.x[-1,0])/Delta; //(\tau_xx)_x
if (y > 0.75)
fprintf(stderr,"Tensor: x %g y %g shear %g normal %g \n",
x, y, shear, normal);
shear = (TVx.y[0,1]+TVx.y[-1,1]-TVx.y[0,-1]-TVx.y[-1,-1])/(4.*Delta);
normal = (TVx.x[1,0]-TVx.x[])/Delta;
if (y > 0.75)
fprintf(stderr,"Vector: x %g y %g shear %g normal %g \n",
x, y, shear, normal);
shear = (TSxy[0,1]+TSxy[-1,1]-TSxy[0,-1]-TSxy[-1,-1])/(4.*Delta);
normal = (TSxx[1,0]-TSxx[])/Delta;
if (y > 0.75)
fprintf(stderr,"Scalar: x %g y %g shear %g normal %g \n",
x, y, shear, normal);
}
foreach()
for (int i = -2; i <= 2; i++)
for (int j = -2; j <= 2; j++)
printf ("%g %g %g %g %g %g %g %g\n", x + i*Delta, y + j*Delta,
S.x.x[i,j], S.y.y[i,j], S.x.y[i,j],
TSxx[i,j], TSyy[i,j], TSxy[i,j]);
}
The results are indeed the same:
Tensor: x 0 y 0.875 shear 4 normal 0
Vector: x 0 y 0.875 shear 4 normal 0
Scalar: x 0 y 0.875 shear 4 normal 0
....