A study of $B_c^+\to J/\psi D_s^+$ and $B_c^+\to J/\psi D_s^{*+}$ decays using 139 fb$^{-1}$ of integrated luminosity collected with the ATLAS detector from $\sqrt{s} = 13$ TeV $pp$ collisions at the LHC is presented. The ratios of the branching fractions of the two decays to the branching fraction of the $B_c^+\to J/\psi \pi^+$ decay are measured: $\mathcal B(B_c^+\to J/\psi D_s^+)/\mathcal B(B_c^+\to J/\psi \pi^+) = 2.76\pm 0.47$ and $\mathcal B(B_c^+\to J/\psi D_s^{*+})/\mathcal B(B_c^+\to J/\psi \pi^+) = 5.33\pm 0.96$. The ratio of the branching fractions of the two decays is found to be $\mathcal B(B_c^+\to J/\psi D_s^{*+})/\mathcal B(B_c^+\to J/\psi D_s^+) = 1.93\pm0.26$. For the $B_c^+\to J/\psi D_s^{*+}$ decay, the transverse polarization fraction, $\Gamma_{\pm\pm}/\Gamma$, is measured to be $0.70\pm0.11$. The reported uncertainties include both the statistical and systematic components added in quadrature. The precision of the measurements exceeds that in all previous studies of these decays. These results supersede those obtained in the earlier ATLAS study of the same decays with $\sqrt{s} = 7$ and 8 TeV $pp$ collision data. A comparison with available theoretical predictions for the measured quantities is presented.