We develop a consistent mathematical
            theory that describes nonlinear interaction of
            wavepackets in a nonlinear periodic dielectric media for
            the dimensions one, two and three. The theory is based
            on the Maxwell equations with quadratic and cubic
            constitutive relations.  Solutions of the Maxwell
            equations on long time intervals are expanded in
            convergent series with respect to a small parameter
            alpha, which measures the contribution of the
            nonlinearity. After that we investigate in detail the
            principal term of the expansion in the case where the
            excitations (and, consequtively, solutions) have a form
            of wavepackets. The ratio of the amplitude frequency and
            the carrier frequency of a wavepacket is an important
            small parameter rho.  The principal term describes the
            nonlinear interaction of a continuum of modes and is
            written in the form of an oscillatory integral. The
            phase function of the integral is written in terms of
            the Floquet-Bloch dispersion relations of the periodic
            media. We consider the situation when the stationary
            phase method is applicable, that means that initially
            the spatial extension of the wavepacket is larger then
            the period and much smaller then the domain where the
            medium is periodic. A detailed mathematical analysis
            shows that the interaction integral expands into sum of
            terms with different powers of rho and the leading terms
            correspond to only a few interacting modes.  We give a
            classification of the nonlinear interactions between
            wavepackets in a media with generic dispersion relations
            based on the powers. The powers take on only a
            relatively small number of prescribed values collected
            in a table, their values depend on a type of degeneracy
            of phase functions formed by the dispersion relations of
            the media. The crucial role in selecting the strongest
            interactions, in particular the second harmonic
            generation in a quadratic medium, is played by internal
            symmetries of the phase function.