/** 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): $$\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 .... ~~~ */