Regularity of harmonic forms on manifolds-with-boundaries

Let $$M$$ be a smooth, bounded, oriented Riemannian manifold-with-boundary. Let $$\alpha$$ be a harmonic differential $$p$$-form on $$M$$, subject to the boundary condition $$\alpha\wedge\nu^\sharp|\partial M = 0$$ or $$\iota_\nu \alpha|\partial M = 0$$. Here $$\nu$$ is the normal vectorfield along $$\partial M$$, $$\nu^\sharp$$ is its dual $$1$$-form, and $$\iota$$ is the interior multiplication. Assume that $$\alpha \in W^{1,2}$$.

The question is: can we conclude that $$\alpha \in C^\infty$$?

The answer depends on your definition of harmonic; if you mean $$\Delta \alpha=0$$ then you can not conclude that $$\alpha$$ is smooth. It is easy to find counter examples on the unit disk in $$\mathbb{R}^2$$. Indeed a $$1$$-form on the unit disk which is a solution of the equation $$\Delta \alpha=0$$ can be written as $$\alpha=f(x,y)dx+g(x,y)dy$$ where $$f$$ and $$g$$ are harmonic functions on the disk. Hence $$f$$ and $$g$$ are uniquely determined by their boundary values on the circle. Let now $$\varphi\colon \partial \mathbb{D}^2\rightarrow \mathbb{R}$$ be a continous but not smooth function and let $$f$$ be the harmonic extension of $$y\varphi$$ and $$g$$ be the harmonic extension of $$-x\varphi$$. Then $$\alpha=f(x,y)dx+g(x,y)dy$$ solves the equation $$\Delta \alpha=0$$ and $$\left.\iota_{\nu} \alpha\right|_{\partial \mathbb{D}^2}=xf+yg=0$$.
If you mean $$(d+d^*) \alpha=0$$ then $$\alpha$$ is smooth because this boundary condition is elliptic (see the chapter 5 in the beautiful book of M. Taylor: Partial Differential Equations I - Basic Theory).