We present a lattice determination of the vector and scalar form factors of the D→πℓν and D→Kℓν semileptonic decays, which are relevant for the extraction of the CKM matrix elements |Vcd| and |Vcs| from experimental data. Our analysis is based on the gauge configurations produced by the European Twisted Mass Collaboration with Nf=2+1+1 flavors of dynamical quarks, at three different values of the lattice spacing (a≃0.062,0.082,0.089  fm) and with pion masses as small as 210 MeV. Quark momenta are injected on the lattice using nonperiodic boundary conditions. The matrix elements of both vector and scalar currents are determined for plenty of kinematical conditions in which parent and child mesons are either moving or at rest. Lorentz symmetry breaking due to hypercubic effects is clearly observed in the data and included in the decomposition of the current matrix elements in terms of additional form factors. After the extrapolations to the physical pion mass and to the continuum limit, we determine the vector and scalar form factors in the whole kinematical region from q2=0 up to qmax2=(MD−Mπ(K))2 accessible in the experiments, obtaining a good overall agreement with experiments, except in the region at high values of q2 where some deviations are visible. A set of synthetic data points, representing our results for f+Dπ(K)(q2) and f0Dπ(K)(q2) for several selected values of q2, is provided and also the corresponding covariance matrix is available. At zero four-momentum transfer, we get f+D→π(0)=0.612(35) and f+D→K(0)=0.765(31). Using the experimental averages for |Vcd|f+D→π(0) and |Vcs|f+D→K(0), we extract |Vcd|=0.2330(137) and |Vcs|=0.945(38), respectively. The second row of the CKM matrix is found to be in agreement with unitarity within the current uncertainties: |Vcd|2+|Vcs|2+|Vcb|2=0.949(78).