diff --git a/ufl/form.py b/ufl/form.py
index 059f17a17..be7d4ae24 100644
--- a/ufl/form.py
+++ b/ufl/form.py
@@ -104,6 +104,12 @@ def coefficients(self):
             self._analyze_form_arguments()
         return self._coefficients
 
+    def geometric_quantities(self):
+        """Return all ``GeometricQuantity`` objects found in form."""
+        if self._geometric_quantities is None:
+            self._analyze_form_arguments()
+        return self._geometric_quantities
+
     def ufl_domain(self):
         """Return the single geometric integration domain occuring in the base form.
 
@@ -595,20 +601,31 @@ def _analyze_domains(self):
         from ufl.geometry import GeometricQuantity
         from ufl.domain import join_domains, sort_domains
 
-        # Collect integration domains
+        # Collect integration domains.
         self._integration_domains = sort_domains(join_domains([itg.ufl_domain() for itg in self._integrals]))
-        # Collect domains in integrands systematically
+        # Collect domains in integrands.
+        """
         domains_in_integrands = []
         for integral in self._integrals:
             integrand = integral.integrand()
             for t in traverse_unique_terminals(integrand):
                 if isinstance(t, (Coefficient, Constant, Argument, GeometricQuantity)):
                     domain = extract_unique_domain(t, expand_mixed_mesh=False)
-                    domains = domain._meshes if isinstance(domain, MixedMesh) else (domain, )
-                    for d in domains:
+                    for d in domain.meshes:
                         if d not in self._integration_domains and d not in domains_in_integrands:
                             domains_in_integrands.append(d)
+        """
+        domains_in_integrands = set()
+        for o in chain(self.arguments(),
+                       self.coefficients(),
+                       self.constants(),
+                       self.geometric_quantities()):
+            domain = extract_unique_domain(o, expand_mixed_mesh=False)
+            domains_in_integrands.update(domain.meshes)
+        domains_in_integrands -= set(self._integration_domains)
         all_domains = self._integration_domains + sort_domains(join_domains(domains_in_integrands))
+        # Let problem solving environments access all domains via
+        # self._domain_numbering.keys() (wrapped in extract_domains()).
         self._domain_numbering = dict((d, i) for i, d in enumerate(all_domains))
 
     def _analyze_subdomain_data(self):