We give a Dirichlet form approach for the construction and analysis of elliptic diffusions in $\bar{\varOmega}\subset\mathbb{R}^n$ with reflecting boundary condition. The problem is formulated in an $L^2$-setting with respect to a reference measure $\mu$ on $\bar{\varOmega}$ having an integrable, $\mathrm{d} x$-almost everywhere (a.e.) positive density $\varrho$ with respect to the Lebesgue measure. The symmetric Dirichlet forms $(\mathcal{E}^{\varrho,a},D(\mathcal{E}^{\varrho,a}))$ we consider are the closure of the symmetric bilinear forms

\begin{gather*} \mathcal{E}^{\varrho,a}(f,g)=\sum_{i,j=1}^n\int_{\varOmega}\partial_ifa_{ij} \partial_jg\,\mathrm{d}\mu,\quad f,g\in\mathcal{D}, \\ \mathcal{D}=\{f\in C(\bar{\varOmega})\mid f\in W^{1,1}_{\mathrm{loc}}(\varOmega),\ \mathcal{E}^{\varrho,a}(f,f)\lt\infty\}, \end{gather*}

in $L^2(\bar{\varOmega},\mu)$, where $a$ is a symmetric, elliptic, $n\times n$-matrix-valued measurable function on $\bar{\varOmega}$. Assuming that $\varOmega$ is an open, relatively compact set with boundary $\partial\varOmega$ of Lebesgue measure zero and that $\varrho$ satisfies the Hamza condition, we can show that $(\mathcal{E}^{\varrho,a},D(\mathcal{E}^{\varrho,a}))$ is a local, quasi-regular Dirichlet form. Hence, it has an associated self-adjoint generator $(L^{\varrho,a},D(L^{\varrho,a}))$ and diffusion process $\bm{M}^{\varrho,a}$ (i.e. an associated strong Markov process with continuous sample paths). Furthermore, since $1\in D(\mathcal{E}^{\varrho,a})$ (due to the Neumann boundary condition) and $\mathcal{E}^{\varrho,a}(1,1)=0$, we obtain a conservative process $\bm{M}^{\varrho,a}$ (i.e. $\bm{M}^{\varrho,a}$ has infinite lifetime). Additionally, assuming that $\sqrt{\varrho}\in W^{1,2}(\varOmega)\cap C(\bar{\varOmega})$ or that $\varrho$ is bounded, $\varOmega$ is convex and $\{\varrho=0\}$ has codimension at least 2, we can show that the set $\{\varrho=0\}$ has $\mathcal{E}^{\varrho,a}$-capacity zero. Therefore, in this case we can even construct an associated conservative diffusion process in $\{\varrho>0\}$. This is essential for our application to continuous $N$-particle systems with singular interactions. Note that for the construction of the self-adjoint generator $(L^{\varrho,a},D(L^{\varrho,a}))$ and the Markov process $\bm{M}^{\varrho,a}$ we do not need to assume any differentiability condition on $\varrho$ and $a$. We obtain the following explicit representation of the generator for $\sqrt{\varrho}\in W^{1,2}(\varOmega)$ and $a\in W^{1,\infty}(\varOmega)$:

$$ L^{\varrho,a}=\sum_{i,j=1}^n\partial_i(a_{ij}\partial_j)+\partial_i(\log\varrho)a_{ij}\partial_j. $$

Note that the drift term can be singular, because we allow $\varrho$ to be zero on a set of Lebesgue measure zero. Our assumptions in this paper even allow a drift that is not integrable with respect to the Lebesgue measure.