This paper proposes a new numerical analysis procedure for reinforced concrete structures with material nonlinearity and cracking. Traditional finite element analysis (FEA), with its reliance on matrix algebra, is an essentially linear procedure. To use it for nonlinear problems, the loading has to be performed in small steps and iterative procedures used to distribute residual nodal forces on each step. The convergence of these procedures degrades after the events of fracturing or cracking in the concrete. If the cracking is extensive, the algorithm may become critically unstable, which may not be equivalent to physical failure. Instead, it is proposed to model the response of reinforced concrete structures in simple tests "constant preload + active load to failure" as a non-step procedure going straight to the target load, To accomplish this, deformation properties of concrete and steel are described by non-associated work-hardening theory of plasticity, formulated such that the plastic strains have a potential. The method is formulated in terms of stresses because the strength criteria are more readily formulated in these terms and because equations of static equilibrium (Navier) are easy to solve. Material nonlinearity of concrete and steel and imperfect bond between them is accounted for by direct method of energy, which is a finite function of stress even if the deformation curve is nonlinear. A variational principle of minimum complementary energy is formulated and a numerical method developed for minimizing it on the domain of statically admissible stress fields (SASF). For construction of SASF for arbitrarily shaped structures in the case of plane stress/plane strain problems, a variant of finite element procedure is proposed. This allows exclusion of "dependent" variables and formulation of the energy functional minimization problem in terms of "fundamental" variables only, which in 2D problems are only one per node. Cracking in the concrete is considered by imposing additional conditions on SASF. For functional minimization, a variant of first-order variable metric method is developed to accommodate "leaps" in the path of descent due to cracking. A methodology for ultimate load capacity analysis is proposed, based on solution of a double optimization (or maximin) problem in the SASF domain. Thus, the failure load is determined as part of the normal computational process and not as a disruption of this process, which improves its credibility.