Functional analysis of variance (ANOVA) modeling has been proved particularly useful to investigate the dynamic changes of functional data according to certain categorical factors and their interactions. However, the current existing methods often encounter difficulties when the functional data are high-dimensional, non-Gaussian, and/or exhibit certain shape characteristics that vary with spatial locations. In this paper, we investigate the functional two-way ANOVA models from a novel Bayesian perspective. A class of highly flexible Gaussian Markov random fields (GMRF) are taken as priors on the functions in the model, which allows us to model various types of functional effects, such as (discrete or continuous) temporal effects and (point-level or areal) spatial effects. The resulting posterior distributions are obtained by an efficient computational tool based on integrated nested Laplace approximations (INLA) (Rue et al., 2009). We then employ the excursion method introduced by Bolin and Lindgren (2015) to build simultaneous credible intervals of functional effects and test their significance from a Bayesian point of view. A simulation study and multiple real data examples are presented to demonstrate the merits of our method.