$\newcommand{\v}[1]{\boldsymbol{#1}}$
First about the appropriate boundary condition we would like to impose: The cross product(scalar-$\mathrm{curl}$) in $\mathbb{R}^2$ is done by we embed a vector into $\mathbb{R}^3$: $\v{v} = \langle v_1,v_2,0\rangle$, then apply the cross product, what we get would be
$$\nabla \times \v{v} = \left\langle 0,0,\left|\begin{matrix}\partial_x&\partial_y\newline v_1&v_2 \end{matrix}\right|\right\rangle
$$
Similarly the vector-$\mathbf{curl}$ is we embed a scalar function into the third dimension of $\mathbb{R}^3$: $u = \langle0,0,u\rangle$, and what we get is the rotation of the $\nabla$-operator normally denoted by $\nabla^{\perp}$:
$$
\nabla^{\perp} u = \left\langle\frac{\partial u}{\partial y},-\frac{\partial u}{\partial x},0 \right\rangle
$$
I read this way of defining 2D-curl operators in a famous computational fluid dynamics book Finite element methods for Navier-Stokes equations: theory and algorithms.
Then in $\mathbb{R}^2$ your problem can be written as:
\begin{equation}
\nabla^{\perp}(\mu^{-1}\nabla \times \v{A}) = \mu \v{J} \tag{1}
\end{equation}
Suppose we multiply a test vector $\v{v}$ vanishing on boundary on both side or the equation and integration by parts:
$$
\int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v}
- \int_{\partial \Omega} \mu^{-1}(\nabla \times \v{A})(\v{v} \times \v{n}) = \int_{\Omega} \mu\v{J}\cdot \v{v}
$$
where recall in $\mathbb{R}^2$, $B = \mu^{-1}(\nabla \times \v{A})$ is really a scalar instead of vector, hence for the tangential Neumann boundary condition, doing the $\mathbb{R}^3$-embedding we have:
$$
\v{n}\times B = \langle n_2 B, -n_1 B,0 \rangle
$$
this is the value of the scalar $B$ along the tangential direction of the domain $\v{t} = \langle n_2, -n_1,0 \rangle$, dot product with $\v{t}$ itself we could get $B$($\v{t}$ is the rotation of a unit outward normal vector), hence the Neumann type boundary condition you really wanna impose here for 2D is:
$$
\mu^{-1}(\nabla \times \v{A}) = g
$$
where $g$ is a scalar function, and your variational form becomes:
$$
\int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v} = \int_{\Omega} \mu\v{J}\cdot \v{v} + \int_{\partial \Omega} g(\v{v} \cdot \v{t}) \tag{2}
$$
The finite element method is we triangulating the domain, construct a finite dimensional Hilbert space $\v{V}_h$ (usually a subspace of the solution space $\v{V}$ to above variational problem), the boundary integral is treated as:
$$
\int_{\partial \Omega} g(\v{v} \cdot \v{t}) = \sum_{e\subset \partial \Omega}
\int_{e} g (\v{v}\cdot \v{t}_e)
$$
where $e$ is a boundary edge, and here is the place quadrature kicks in, that you wanna evaluate a line integral on each edge, it depends on what your finite element test function space $\v{V}_h$, nowadays a consensus is to choose Nédélec elements, in lowest order case it is $\v{v}_e = \lambda_A\nabla \lambda_B - \lambda_B\nabla \lambda_A $ for $e = \overrightarrow{AB}$.
Secondly, what is really problematic with your statement here is:
for a smooth vector field $\v{A} \in \v{C}^{\infty}(\mathbb{R}^2)$ or $\v{C}^{\infty}_c(\Omega)$, the equation (1) doesn't simply reduce to the scalar Laplace equation as you mentioned in OP, a direct observation would be that homogeneous Neumann problem for $-\Delta u = f$
has a constant as null space because a solution plus a constant is still a solution, while $\v{curl}$-problem has all the gradient vector fields as null space!
Similar thing applies for variational formulation (2), the bilinear form is not coercive, simply plugging a gradient field $\v{A} = \v{v} = \nabla \varphi$ into it. What we numerical analyst usually do here is to find a solution in the quotient space $\v{H}(\mathrm{curl})/\nabla H^1$ by adding a Lagrange multiplier to let the problem fall into mixed finite element(google MFEM if you would like to know more) framework:
$$
\begin{cases}
\int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v}
+ \int_{\Omega} \v{v}\cdot \nabla \varphi= \int_{\Omega} \mu\v{J}\cdot \v{v} + \int_{\partial \Omega} g(\v{v} \cdot \v{t}) & \text{ for } \forall \v{v}\in \v{H}(\mathrm{curl})\newline
\int_{\Omega} \v{A} \cdot \nabla \psi = 0& \text{ for } \forall \psi\in H^1/\mathbb{R}
\end{cases}
$$
to make sure problem (2) is solvable using finite element, not only we have to use Lagrange multiplier to circumvent the "null space" problem, we also have to guarantee that the right side is in the range of the curl-operator, here plug $\v{v} = \nabla \varphi$ into (2):
$$
0=\int_{\Omega} \mu\v{J}\cdot \nabla \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \v{t}) = -\int_{\Omega} \nabla \cdot(\mu\v{J}) \varphi + \int_{\partial \Omega} (\mu\v{J} \cdot \v{n}) \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \v{t})
$$
in your case, $g=0$, also by the arbitrariness of $\varphi$, for eg we choose $\nabla \varphi \cdot \v{t} = 0$ for $\varphi\in H^1_0$, we would get so called "compatibility condition" in finite element method for Neumann problem:
$$
\nabla \cdot(\mu\v{J}) = 0\text{ in }\Omega\quad \text{ and }\quad \mu\v{J}\cdot \v{n} = 0\text{ on }\partial\Omega
$$
Lastly about the implementation of finite element method for this equation, if you have a pure Neumann problem with homogeneous boundary condition(derivative on boundary is zero), then you do not have to do anything during assembling of the stiffness matrix or right hand side integral.
I think you are asking about the possibility of satisfying the desired boundary conditions by each solution of the Helmholtz equation in the product form $X(x) Y(y)$, where $X(x)$ an $Y(y)$ were obtained by separation of variables. This can only happen if the mode functions $X(x)$ and $Y(y)$ have zero sets that are compatible with the chosen boundary. Or more specifically, if the coordinates $x$ and $y$ are compatible with a chosen boundary. For instance, using separation of variables in Cartesian coordinates to solve a boundary value problem with a circular boundary is hopeless, because neither of the $x$ or $y$ coordinates is constant even along any piece of the boundary. So, as you remarked, to use this method, one needs to choose a coordinate system that is compatible both with the differential equation and with the boundary.
You can see the complete list of 11 coordinate systems in which the Helmholtz equation is separable here on MathWorld. I think your example $|y|=ax$ with $a\ne \pm 1$ is not appropriate for any of these coordinate system. But the $y=ax^2$ example may be appropriate for parabolic coordinates.
Best Answer
If you don't mind wasting some resources, extend the $y$ coordinate to the interval $[0,4\pi]$. Then only look at solutions that are odd (anti-symmetric) under the flip $y\mapsto 4\pi-y$. In particular, if your PDE has source terms, they have to be anti-symmetrically mirrored from the $[0,2\pi]$ to the $[2\pi,4\pi]$ half of the extended interval. When restricted to $y\in[0,2\pi]$ a periodic solution on the extended interval will satisfy the desired boundary conditions.