Local modifications of a computational domain are often performed in order to simplify the meshing process and to reduce computational costs and memory requirements. However, removing geometrical features of a domain often introduces a non-negligible error in the solution of a differential problem in which it is defined. In this paper, we aim at generalizing the work from [1], in which an a posteriori estimator of the geometrical defeaturing error is derived for domains from which one geometrical feature is removed. More precisely, we study the case of domains containing an arbitrary number of distinct features, and we perform an analysis on Poisson's, linear elasticity, and Stokes' equations. We introduce a simple and computationally cheap a posteriori estimator of the geometrical defeaturing error, whose reliability and efficiency are rigorously proved, and we introduce a geometric refinement strategy that accounts for the defeaturing error: Starting from a fully defeatured geometry, the algorithm determines at each iteration step which features need to be added to the geometrical model to reduce the defeaturing error. These important features are then added to the (partially) defeatured geometrical model at the next iteration, until the solution attains a prescribed accuracy. A wide range of numerical experiments are finally reported to illustrate and validate this work.