Let $\Phi_t(x)$ be the flow of your vector field $\nabla d_{\mathcal{A}}$, starting at $x \in \mathcal{A}$. It will exist until the time $t_0(x)$, when it hits the set where $\nabla d_{\mathcal{A}}$ is not smooth.
Consider the ODE
$$2 \frac{\partial}{\partial t} z(t, x) + z(t, x)\cdot \Delta d_{\mathcal{A}}(\Phi_t(x)) = 0, ~~~~~~~ z(0, x) = W(x)$$
depending on the parameter $x \in \partial \mathcal{A}$. It is a linear ODE, hence it has a unique solution that exists for all time, or rather until $t_0(x)$ (because afterwards the equation does not make sense anymore). It is also well-known that the solution depends smoothly on the parameter $x$, as all the data is smooth.
Now let $\mathcal{B} = \{\Phi_t(x) \mid x \in \partial \mathcal{A}, t < t_0(x)\} \subseteq \mathcal{A}$ be the set of points $y$ that lie on a unique flowline, i.e. there exist a unique $t$ and a unique $x \in \partial \mathcal{A}$ such that $y = \Phi_t(x)$ (that this set has the claimed description is because flowlines cannot cross). This is a neighborhood of $\partial \mathcal{A}$ in $\mathcal{A}$, and its complement in $\mathcal{A}$ has measure zero.
For points in $\mathcal{B}$, set
$$Z(\Phi_t(x)) := z(t, x).$$
This gives you a solution that has the desired boundary values and is smooth on $\mathcal{B}$.
However, you know nothing about the smoothness in $\mathcal{A} \setminus \mathcal{B}$. For example, if you take a disc in $\mathbb{R}^n$, your function $Z$ will be continuous if and only if $W(x) \equiv \mathrm{const}$. Otherwise it will not be continuous at the origin.
The Neumann boundary condition is a "natural" BC. You don't need to impose it. You have to change the space $H^1_0$ for $H^1$ and you end up with the same weak formulation that you found in the Dirichlet BC case, only that you required that the test functions lie in $H^1$.
That will do it.
Best Answer
A weak form of your BVP is $a(q,v)=\ell(v)$ where $a(q,v)=\int_{\Omega}\nabla q\cdot\nabla v\,dx+\int_{\Omega}\beta qv\,dx$ and $\ell(v)=\int_{\Omega}fv\,dx+\int_{\partial\Omega}gv\,ds$ with $q,v\in H^1(\Omega)$. If $\beta>0$, the bilinear form is coercive and continuous in $H^1(\Omega)$. Thus, apply Lax-Milgram and you get existence, uniqueness and stability in $H^1(\Omega)$. Stability now depends on the value of $\beta$.