We consider the nonlinear problem of anisotropic Allen-Cahn equation \begin{document}$ \varepsilon^2{\mathrm {div}}\big( \nabla_{{\mathfrak a}(y)} u\big)+\mathcal{P}(y)u(1-u^2) = 0\quad \mbox{in}\ \Omega, \qquad \nabla_{{\mathfrak a}(y)} u\cdot \nu = 0\quad \mbox{on}\ \partial \Omega, $\end{document} where \begin{document}$ \Omega $\end{document} is a bounded domain in \begin{document}$ {\mathbb R}^2 $\end{document} with smooth boundary, \begin{document}$ \varepsilon $\end{document} is a small positive parameter, \begin{document}$ \nu $\end{document} denotes the unit outward normal of \begin{document}$ \partial\Omega $\end{document}, and \begin{document}$ \mathcal{P}(y) $\end{document} is a uniformly positive smooth potential on \begin{document}$ \bar{\Omega} $\end{document}. The operator \begin{document}$ \nabla_{{\mathfrak a}(y)} u $\end{document} is defined by$ \nabla_{{\mathfrak a}(y)} u = \big({\mathfrak a}_1(y)u_{y_1}, {\mathfrak a}_2(y)u_{y_2}\big) \quad \mbox{with }\ {\mathfrak a}(y) = \big({\mathfrak a}_1(y), {\mathfrak a}_2(y)\big), $where \begin{document}$ {\mathfrak a}_1(y) $\end{document} and \begin{document}$ {\mathfrak a}_2(y) $\end{document} are two positive smooth functions on \begin{document}$ \bar\Omega $\end{document}. Let \begin{document}$ \Gamma $\end{document} be an interior curve intersecting orthogonally with \begin{document}$ \partial\Omega $\end{document} at exactly two points or a closed simple curve in \begin{document}$ \Omega $\end{document}, and dividing \begin{document}$ \Omega $\end{document} into two parts. Moreover, \begin{document}$ \Gamma $\end{document} is a non-degenerate geodesic embedded in the Riemannian manifold \begin{document}$ {\mathbb R}^2 $\end{document} associated with metric \begin{document}$ \mathcal{P}(y)\big({\mathfrak a}_2(y){\mathrm d}y_1\otimes{\mathrm d}y_1+{\mathfrak a}_1(y){\mathrm d}y_2\otimes{\mathrm d}y_2\big) $\end{document}. By assuming some additional constraints on the functions \begin{document}$ {\mathfrak a}(y) $\end{document}, \begin{document}$ \mathcal{P}(y) $\end{document} and the curves \begin{document}$ \Gamma $\end{document}, \begin{document}$ \partial\Omega $\end{document}, we prove that there exists a solution \begin{document}$ u_{\varepsilon} $\end{document} with an interface such that: as \begin{document}$ \varepsilon\rightarrow 0 $\end{document}, \begin{document}$ u_{\varepsilon} $\end{document} approaches \begin{document}$ +1 $\end{document} in one part of \begin{document}$ \Omega $\end{document}, while tends to \begin{document}$ -1 $\end{document} in the other part, except a small neighborhood of \begin{document}$ \Gamma $\end{document}.