We present a novel approach for analysing multivariate case‐control georeferenced data in a Bayesian disease mapping context using stochastic partial differential equations (SPDEs) and the integrated nested Laplace approximation (INLA) for model fitting. In particular, we propose smooth terms based on SPDE models to estimate the underlying spatial variation as well as risk associated to pollution sources. Log‐Gaussian Cox processes are used to estimate the intensity of the cases and controls, to account for risk factors and include a term to measure spatial residual variation. Each intensity is modelled on a baseline spatial effect (estimated from both controls and cases), a disease‐specific spatial term and the effects of some covariates. By fitting these models, the residual spatial terms can be easily compared to detect high‐risk areas not explained by the covariates. Three different types of effects to model exposure to pollution sources are considered on the distance to the source: a fixed effect, a smooth term to model non‐linear effects by means of a discrete random walk of order one and a Gaussian process in one dimension with a Matérn covariance function. Spatial terms are modelled using a Gaussian process in two dimensions with a Matérn covariance function and are approximated using an approach based on solving an SPDE through INLA. Finally, this new framework is applied to a dataset of three different types of cancer and a set of controls from Alcalá de Henares (Madrid, Spain). Covariates available include the distance to several polluting industries and socioeconomic indicators. Our findings point to a possible risk increase due to the proximity to some of these industries.