The production of dark matter in association with Higgs bosons is predicted in several extensions of the Standard Model. An exploration of such scenarios is presented, considering final states with missing transverse momentum and $b$-tagged jets consistent with a Higgs boson. The analysis uses proton-proton collision data at a centre-of-mass energy of 13 TeV recorded by the ATLAS experiment at the LHC during Run 2, amounting to an integrated luminosity of 139 fb$\mathrm{<mrow\; data-mjx-texclass="ORD">-1</mrow>}$. The analysis, when compared with previous searches, benefits from a larger dataset, but also has further improvements providing sensitivity to a wider spectrum of signal scenarios. These improvements include both an optimised event selection and advances in the object identification, such as the use of the likelihood-based significance of the missing transverse momentum and variable-radius track-jets. No significant deviation from Standard Model expectations is observed. Limits are set, at 95% confidence level, in two benchmark models with two Higgs doublets extended by either a heavy vector boson $Z\prime$ or a pseudoscalar singlet $a$ and which both provide a dark matter candidate $\chi$. In the case of the two-Higgs-doublet model with an additional vector boson $Z\prime$, the observed limits extend up to a $Z\prime$ mass of 3 TeV for a mass of 100 GeV for the dark matter candidate. The two-Higgs-doublet model with a dark matter particle mass of 10 GeV and an additional pseudoscalar $a$ is excluded for masses of the $a$ up to 520 GeV and 240 GeV for $\tan \beta =1$ and $\tan \beta =10$ respectively. Limits on the visible cross-sections are set and range from 0.05 fb to 3.26 fb, depending on the missing transverse momentum and $b$-quark jet multiplicity requirements.