Solving a Stress Problem with VECFEM
The vector of distortions in global cartesian coordinates is defined by

       E1  = U1x
       E2  = U2y
       E3  = U3z
       E12 = (U2x+U1y)/2
       E23 = (U3y+U2z)/2
       E13 = (U3x+U1z)/2.

The corresponding stress vector has the form (S1,S2,S3,S12,S23,S13). By
Hooke's law,  the stress and distortion vector for orthotropic
materials corresponds as follows

    | S1  |       | C11  C12 C13  0     0    0   |   | E1  |
    | S2  |       | C12  C22 C23  0     0    0   |   | E2  |
    | S3  |   =   | C13  C23 C33  0     0    0   | * | E3  |
    | S12 |       |  0    0   0 2*C44   0    0   |   | E12 |
    | S23 |       |  0    0   0   0   2*C55  0   |   | E23 |
    | S13 |       |  0    0   0   0     0  2*C66 |   | E13 |

The matrix is called the tensor of elasticity. For the case of isotropic
material the C-enteries are given by the two values E and NY called
modulus of elasticity and Passion's number.

The inner energy of the body tracted by the nodal loades
(F1(i),F2(i),F3(i)) at the nodes X(i)=(x(i),y(i),z(i)) (i=1,NFOR) and
distorted by the displacements U is given by

       1/2 * L(U,U)-F(U)
where L is a certain bilinear form and F is a linear form. The bilinear
form L is an integration over the body and consists of the interactions
of the derivative of the test functions and the derivative of the
solution. The following table gives the entries in the L3-matrix
(L2=L1=L0=0). The rows are the test function (the left number counts the
component of the test function (COMP1) and the right number counts the
space direction of the derivative) and the columns are the solution (the
upper number counts the component of the solution (COMP2) and the lower
number counts the space direction of the derivative) :

      L3
       \  COMP2 I      1      I      2      I      3      I
      COMP1  \  I  1   2   3  I  1   2   3  I  1   2   3  I
      -----------------------------------------------------
             1  I C11  0   0  I  0  C12  0  I  0   0  C13 I
         1   2  I  0  C66  0  I C66  0   0  I  0   0   0  I
             3  I  0   0  C55 I 0   0   0   I C55  0   0  I
      -----------------------------------------------------
             1  I  0  C66  0  I C66  0   0  I  0   0   0  I
         2   2  I C12  0   0  I  0  C22  0  I  0   0  C23 I
             3  I  0   0   0  I  0   0  C44 I  0  C44  0  I
      -----------------------------------------------------
             1  I  0   0  C55 I 0   0   0   I C55  0   0  I
         2   2  I  0   0   0  I  0   0  C44 I  0  C44  0  I
             3  I C13  0   0  I  0  C23  0  I  0   0  C33 I
      -----------------------------------------------------

L is symmetric. The right hand side F is the sum of the evaluations of
the test function at the working points of the concentrated loads
weighted by the value of the forces  (i=1,NFOR) :

      V1(X(i)) * F1(i) + V2(X(i)) * F2(i) + V3(X(i)) * F3(i).

So these points are a 0-dimensional manifold. The values of F0 for the
i-th node in this 0-dimensional manifold is given by the following
table :

        COMP I  F0
        ------------
          1  I F1(i)
          2  I F2(i)
          3  I F3(i)

where the row counts the component of the test function.
The displacements in certain points of the body are prescribed.
Mathematically this is described by Dirichlet conditions.

By the equilibrium condition the real displacement U which fulfills the
prescribed displacements is given by the functional equation

              L(V,U)=F(V)

for all test function V which satify the homogenous Dirichlet
conditions. This functional equation can be solved by VECFEM.