The mathematical modeling by FEM of the piezoelectric devices is considered. By applying semi-discrete FEM approximations of the solution to the governing equations in weak form the variational FEM equations with symmetrical saddle matrices are derived. On the base of symmetrical algorithms the computer program ACELAN were developed. For one- and two-dimensional piezoelectric structures new quasielastic finite elements are obtained. New models for describing process of polarization or repolarization of piezoceramic are suggested. The 3D mathematical models of porous and polycrystalline piezocomposites with 3-0 and 3-3 connectivity are built. The finite element programs for the effective constants and effective properties calculations are developed.