SubDomainSet broken (both functionally and performance)
Consider the following MFE:

```
        class DummySubdomains(SubDomainSet):
            name = 'dummydomain'
        dummy = DummySubdomains(N=1, bounds=(1, 1, 1, 1))

        grid = Grid(shape=(10, 10), subdomains=(dummy,))

        f = TimeFunction(name='f', grid=grid)
        g = TimeFunction(name='g', grid=grid)
        theta = Function(name='theta', grid=grid)
        phi = Function(name='phi', grid=grid)

        eqns = [Eq(f.forward, f*sin(phi), subdomain=grid.subdomains['dummydomain']),
                Eq(g.forward, g*sin(theta), subdomain=grid.subdomains['dummydomain'])]

        op = Operator(eqns)
```

There are three issues:

1) functional: as @mloubout pointed out on Slack , this generates invalid code as the thickness assignments (eg `xi_n_ltkn = xi_n_m[n]`) are generated N times (with N number of equations, here N=2) . Not only are these assignments redundant, but also flawed since the first one assigns to a `const int` variable, which is later assigned to again, thus causing a compilation failure
2) the generation of redundant assignments prevents the fusion of the two eqns into the same loop (performance issue)
3) the hoisted temporaries are redundantly computed within the time loop, while they should go outside and be computed over the whole grid once and for all 
diff --git a/devito/core/cpu.py b/devito/core/cpu.py
index ee1e5eb1c..89df6fa43 100644
--- a/devito/core/cpu.py
+++ b/devito/core/cpu.py
@@ -5,7 +5,7 @@ from devito.exceptions import InvalidOperator
 from devito.passes.equations import collect_derivatives
 from devito.passes.clusters import (Lift, blocking, buffering, cire, cse,
                                     extract_increments, factorize, fission, fuse,
-                                    optimize_pows)
+                                    optimize_pows, optimize_msds)
 from devito.passes.iet import (CTarget, OmpTarget, avoid_denormals, linearize, mpiize,
                                optimize_halospots, hoist_prodders, relax_incr_dimensions)
 from devito.tools import timed_pass
@@ -161,6 +161,9 @@ class Cpu64AdvOperator(Cpu64OperatorMixin, CoreOperator):
         platform = kwargs['platform']
         sregistry = kwargs['sregistry']
 
+        # Optimize MultiSubDomains
+        clusters = optimize_msds(clusters)
+
         # Toposort+Fusion (the former to expose more fusion opportunities)
         clusters = fuse(clusters, toposort=True)
 
@@ -247,6 +250,9 @@ class Cpu64FsgOperator(Cpu64AdvOperator):
         platform = kwargs['platform']
         sregistry = kwargs['sregistry']
 
+        # Optimize MultiSubDomains
+        clusters = optimize_msds(clusters)
+
         # Toposort+Fusion (the former to expose more fusion opportunities)
         clusters = fuse(clusters, toposort=True)
 
diff --git a/devito/core/gpu.py b/devito/core/gpu.py
index 7b76feb7d..ace248513 100644
--- a/devito/core/gpu.py
+++ b/devito/core/gpu.py
@@ -7,7 +7,7 @@ from devito.exceptions import InvalidOperator
 from devito.passes.equations import collect_derivatives
 from devito.passes.clusters import (Lift, Streaming, Tasker, blocking, buffering,
                                     cire, cse, extract_increments, factorize,
-                                    fission, fuse, optimize_pows)
+                                    fission, fuse, optimize_pows, optimize_msds)
 from devito.passes.iet import (DeviceOmpTarget, DeviceAccTarget, optimize_halospots,
                                mpiize, hoist_prodders, is_on_device, linearize)
 from devito.tools import as_tuple, timed_pass
@@ -150,6 +150,9 @@ class DeviceAdvOperator(DeviceOperatorMixin, CoreOperator):
         platform = kwargs['platform']
         sregistry = kwargs['sregistry']
 
+        # Optimize MultiSubDomains
+        clusters = optimize_msds(clusters)
+
         # Toposort+Fusion (the former to expose more fusion opportunities)
         clusters = fuse(clusters, toposort=True, options=options)
 
diff --git a/devito/ir/clusters/algorithms.py b/devito/ir/clusters/algorithms.py
index 2db0d7d2a..819ad9013 100644
--- a/devito/ir/clusters/algorithms.py
+++ b/devito/ir/clusters/algorithms.py
@@ -286,8 +286,10 @@ class Stepper(Queue):
                 exprs = [e.apply(func) for e in exprs]
 
             # Augment IterationSpace
-            ispace = IterationSpace(c.ispace.intervals,
-                                    {**c.ispace.sub_iterators, **{d: tuple(mds)}},
+            sub_iterators = dict(c.ispace.sub_iterators)
+            sub_iterators[d] = tuple(i for i in sub_iterators[d] + tuple(mds)
+                                     if i not in subiters)
+            ispace = IterationSpace(c.ispace.intervals, sub_iterators,
                                     c.ispace.directions)
 
             processed.append(c.rebuild(exprs=exprs, ispace=ispace))
diff --git a/devito/ir/equations/algorithms.py b/devito/ir/equations/algorithms.py
index a0348edaf..072168f66 100644
--- a/devito/ir/equations/algorithms.py
+++ b/devito/ir/equations/algorithms.py
@@ -82,7 +82,7 @@ def generate_implicit_exprs(expressions):
 
     Currently, implicit expressions stem from the following:
 
-        * ``SubDomainSet``'s attached to input equations.
+        * MultiSubDomains attached to input equations.
     """
     found = {}
     processed = []
@@ -91,9 +91,9 @@ def generate_implicit_exprs(expressions):
             try:
                 dims = [d.root for d in e.free_symbols if isinstance(d, Dimension)]
                 sub_dims = [d.root for d in e.subdomain.dimensions]
-                sub_dims.append(e.subdomain.implicit_dimension)
+                sub_dims.extend(e.subdomain.implicit_dimensions)
                 dims = [d for d in dims if d not in frozenset(sub_dims)]
-                dims.append(e.subdomain.implicit_dimension)
+                dims.extend(e.subdomain.implicit_dimensions)
                 if e.subdomain not in found:
                     grid = list(retrieve_functions(e, mode='unique'))[0].grid
                     found[e.subdomain] = [i.func(*i.args, implicit_dims=dims) for i in
@@ -103,6 +103,7 @@ def generate_implicit_exprs(expressions):
                 new_e = Eq(e.lhs, e.rhs, subdomain=e.subdomain, implicit_dims=dims)
                 processed.append(new_e)
             except AttributeError:
+                # Not a MultiSubDomain
                 processed.append(e)
         else:
             processed.append(e)
diff --git a/devito/ir/support/space.py b/devito/ir/support/space.py
index a8aef97b6..f1ca5ce0a 100644
--- a/devito/ir/support/space.py
+++ b/devito/ir/support/space.py
@@ -815,6 +815,15 @@ class IterationSpace(Space):
 
         return IterationSpace(intervals, self.sub_iterators, self.directions)
 
+    def relaxed(self, dims):
+        f = lambda d: d in as_tuple(dims)
+
+        intervals = [i.relaxed if f(i.dim) else i for i in self.intervals]
+        sub_iterators = {k.root if f(k) else k: v for k, v in self.sub_iterators.items()}
+        directions = {k.root if f(k) else k: v for k, v in self.directions.items()}
+
+        return IterationSpace(intervals, sub_iterators, directions)
+
     def is_compatible(self, other):
         """
         A relaxed version of ``__eq__``, in which only non-derived dimensions
diff --git a/devito/ir/support/utils.py b/devito/ir/support/utils.py
index 564c44602..cd9bb6cb6 100644
--- a/devito/ir/support/utils.py
+++ b/devito/ir/support/utils.py
@@ -78,7 +78,7 @@ def build_iterators(mapper):
     iterators = OrderedDict()
     for k, v in mapper.items():
         for d in v:
-            if d.is_Stepping:
+            if d.is_Stepping or d.is_Incr:
                 values = iterators.setdefault(d.root, [])
                 if d not in values:
                     values.append(d)
@@ -98,7 +98,7 @@ def build_intervals(stencil):
     for d, offs in stencil.items():
         if d.is_Stepping:
             mapper[d.root].update(offs)
-        elif d.is_Conditional:
+        elif d.is_Conditional or d.is_Incr:
             mapper[d.parent].update(offs)
         elif d.is_Modulo:
             mapper[d.root].update({d.offset - d.root + i for i in offs})
diff --git a/devito/passes/clusters/aliases.py b/devito/passes/clusters/aliases.py
index 11e58a048..d43b15c88 100644
--- a/devito/passes/clusters/aliases.py
+++ b/devito/passes/clusters/aliases.py
@@ -18,6 +18,7 @@ from devito.tools import (Stamp, as_mapper, as_tuple, flatten, frozendict, gener
                           split)
 from devito.types import (Array, TempFunction, Eq, Symbol, ModuloDimension,
                           CustomDimension, IncrDimension, Indexed)
+from devito.types.grid import MultiSubDimension
 
 __all__ = ['cire']
 
@@ -108,6 +109,7 @@ class CireTransformer(object):
         self.opt_rotate = options['cire-rotate']
         self.opt_ftemps = options['cire-ftemps']
         self.opt_mingain = options['cire-mingain']
+        self.opt_multisubdomain = True
 
     def _aliases_from_clusters(self, clusters, exclude, meta):
         exprs = flatten([c.exprs for c in clusters])
@@ -140,6 +142,10 @@ class CireTransformer(object):
         # Schedule -> [Clusters]_k
         processed, subs = lower_schedule(schedule, meta, self.sregistry, self.opt_ftemps)
 
+        # [Clusters]_k -> [Clusters]_k (optimization)
+        if self.opt_multisubdomain:
+            processed = optimize_clusters_msds(processed)
+
         # [Clusters]_k -> [Clusters]_{k+n}
         for c in clusters:
             n = len(c.exprs)
@@ -315,6 +321,7 @@ class CireSops(CireTransformer):
 
         self.opt_maxpar = options['cire-maxpar']
         self.opt_schedule_strategy = options['cire-schedule']
+        self.opt_multisubdomain = False
 
     def process(self, clusters):
         processed = []
@@ -898,6 +905,40 @@ def lower_schedule(schedule, meta, sregistry, ftemps):
     return clusters, subs
 
 
+def optimize_clusters_msds(clusters):
+    """
+    Relax the clusters by letting the expressions defined over MultiSubDomains to
+    rather be computed over the entire domain. This increases the likelihood of
+    code lifting by later passes.
+    """
+    processed = []
+    for c in clusters:
+        msds = [d for d in c.ispace.itdimensions if isinstance(d, MultiSubDimension)]
+
+        if msds:
+            mapper = {d: d.root for d in msds}
+            exprs = [uxreplace(e, mapper) for e in c.exprs]
+
+            ispace = c.ispace.relaxed(msds)
+
+            accesses = detect_accesses(exprs)
+            parts = {k: IntervalGroup(build_intervals(v)).relaxed
+                     for k, v in accesses.items() if k}
+            intervals = [i for i in c.dspace if i.dim not in msds]
+            dspace = DataSpace(intervals, parts)
+
+            guards = {mapper.get(d, d): v for d, v in c.guards.items()}
+            properties = {mapper.get(d, d): v for d, v in c.properties.items()}
+            syncs = {mapper.get(d, d): v for d, v in c.syncs.items()}
+
+            processed.append(c.rebuild(exprs=exprs, ispace=ispace, dspace=dspace,
+                                       guards=guards, properties=properties, syncs=syncs))
+        else:
+            processed.append(c)
+
+    return processed
+
+
 def pick_best(variants, schedule_strategy, eval_variants_delta):
     """
     Return the variant with the best trade-off between operation count
diff --git a/devito/passes/clusters/misc.py b/devito/passes/clusters/misc.py
index 7dd0863d8..f65b0f7d9 100644
--- a/devito/passes/clusters/misc.py
+++ b/devito/passes/clusters/misc.py
@@ -7,9 +7,10 @@ from devito.passes.clusters.utils import cluster_pass
 from devito.symbolics import pow_to_mul
 from devito.tools import DAG, Stamp, as_tuple, flatten, frozendict, timed_pass
 from devito.types import Symbol
+from devito.types.grid import MultiSubDimension
 
 __all__ = ['Lift', 'fuse', 'optimize_pows', 'extract_increments',
-           'fission']
+           'fission', 'optimize_msds']
 
 
 class Lift(Queue):
@@ -365,3 +366,68 @@ def fission(clusters):
                                    ..
     """
     return Fission().process(clusters)
+
+
+class MSDOptimizer(Queue):
+
+    """
+    Implement MultiSubDomains optimization.
+
+    Currently, the following optimizations are performed:
+
+        * Removal of redundant thicknesses assignments. These stem from Eqs
+          defined over the same MultiSubDomain in the very same loop nest.
+          The redundant assignments obviously do not impact correctness,
+          but they may affect other optimizations, such as fusion.
+    """
+
+    def callback(self, clusters, prefix):
+        if not prefix or any(isinstance(i.dim, MultiSubDimension) for i in prefix):
+            return clusters
+
+        msds = {d for d in set().union(*[c.dimensions for c in clusters])
+                if isinstance(d, MultiSubDimension)}
+        if not msds:
+            return clusters
+
+        # Remove redundant thicknesses assignments
+
+        thicknesses = set().union(*[list(i._thickness_map) for i in msds])
+        candidates = [c for c in clusters if set(c.scope.writes) & thicknesses]
+
+        # First of all, make sure we analyze all and only the thicknesses assignments
+        # at the same depth
+        d = prefix[-1].dim
+        if any(c.itintervals[-1].dim is not d for c in candidates):
+            return clusters
+
+        # Then, attempt extirpation of redundancies
+        schedulable = set(thicknesses)
+        processed = []
+        for c in clusters:
+            if c in candidates:
+                exprs = []
+                for e in c.exprs:
+                    try:
+                        schedulable.remove(e.lhs)
+                        exprs.append(e)
+                    except KeyError:
+                        # Already scheduled, no-op
+                        pass
+                if exprs:
+                    processed.append(c.rebuild(exprs=exprs))
+            else:
+                processed.append(c)
+
+        # Sanity check
+        assert len(schedulable) == 0
+
+        return processed
+
+
+@timed_pass()
+def optimize_msds(clusters):
+    """
+    Optimize clusters defined over MultiSubDomains.
+    """
+    return MSDOptimizer().process(clusters)
diff --git a/devito/types/grid.py b/devito/types/grid.py
index 9731daac5..821b3e534 100644
--- a/devito/types/grid.py
+++ b/devito/types/grid.py
@@ -461,6 +461,10 @@ class SubDomain(object):
     def dimensions(self):
         return self._dimensions
 
+    @property
+    def ndim(self):
+        return len(self.dimensions)
+
     @property
     def dimension_map(self):
         return {d.root: d for d in self.dimensions}
@@ -480,6 +484,24 @@ class SubDomain(object):
         """
         raise NotImplementedError
 
+    def union(self, other):
+        """
+        Return the union of two SubDomains as a new SubDomain.
+        """
+        raise NotImplementedError
+
+    def intersection(self, other):
+        """
+        Return the intersection of two SubDomains as a new SubDomain.
+        """
+        raise NotImplementedError
+
+    def difference(self, other):
+        """
+        Return the difference of two SubDomains as a new SubDomain.
+        """
+        raise NotImplementedError
+
 
 class Domain(SubDomain):
 
@@ -505,7 +527,88 @@ class Interior(SubDomain):
         return {d: ('middle', 1, 1) for d in dimensions}
 
 
-class SubDomainSet(SubDomain):
+class MultiSubDimension(SubDimension):
+
+    """
+    A special SubDimension for graceful lowering of MultiSubDomains.
+    """
+
+    def __init_finalize__(self, name, parent):
+        lst, rst = self._symbolic_thickness(name)
+        left = parent.symbolic_min + lst
+        right = parent.symbolic_max - rst
+
+        super().__init_finalize__(name, parent, left, right, ((lst, 0), (rst, 0)), False)
+
+
+class MultiSubDomain(SubDomain):
+
+    """
+    Base class for types representing a group of SubDomains.
+    """
+
+    @classmethod
+    def _bounds_glb_to_loc(cls, dec, m, M):
+        """
+        Translate a SubDomain global bounds, that is thicknesses, into local bounds.
+        """
+        bounds_m = np.zeros(m.shape, dtype=m.dtype)
+        bounds_M = np.zeros(m.shape, dtype=m.dtype)
+        for j in range(m.size):
+            lmin = dec.glb_min + m[j]
+            lmax = dec.glb_max - M[j]
+
+            # Check if the subdomain doesn't intersect with the decomposition
+            if lmin < dec.loc_abs_min and lmax < dec.loc_abs_min:
+                bounds_m[j] = dec.loc_abs_max
+                bounds_M[j] = dec.loc_abs_max
+                continue
+            if lmin > dec.loc_abs_max and lmax > dec.loc_abs_max:
+                bounds_m[j] = dec.loc_abs_max
+                bounds_M[j] = dec.loc_abs_max
+                continue
+
+            if lmin < dec.loc_abs_min:
+                bounds_m[j] = 0
+            elif lmin > dec.loc_abs_max:
+                bounds_m[j] = dec.loc_abs_max
+                bounds_M[j] = dec.loc_abs_max
+                continue
+            else:
+                bounds_m[j] = dec.index_glb_to_loc(m[j], LEFT)
+
+            if lmax < dec.loc_abs_min:
+                bounds_m[j] = dec.loc_abs_max
+                bounds_M[j] = dec.loc_abs_max
+                continue
+            elif lmax >= dec.loc_abs_max:
+                bounds_M[j] = 0
+            else:
+                bounds_M[j] = dec.index_glb_to_loc(M[j], RIGHT)
+
+        return bounds_m, bounds_M
+
+    def __subdomain_finalize__(self, dimensions, shape, distributor=None, **kwargs):
+        """
+        Must be overridded by subclasses.
+        """
+        raise NotImplementedError
+
+    def _create_implicit_exprs(self, grid):
+        """
+        Must be overridded by subclasses.
+        """
+        raise NotImplementedError
+
+    @property
+    def implicit_dimensions(self):
+        """
+        The Dimensions implicitly defined by the MultiSubDomain.
+        """
+        return ()
+
+
+class SubDomainSet(MultiSubDomain):
     """
     Class to define a set of N (a positive integer) subdomains.
 
@@ -583,22 +686,20 @@ class SubDomainSet(SubDomain):
     """The implicit dimension of the SubDomainSet."""
 
     def __init__(self, **kwargs):
-        super(SubDomainSet, self).__init__()
+        super().__init__()
+
         if self.implicit_dimension is None:
-            n = Dimension(name='n')
-            self.implicit_dimension = n
+            self.implicit_dimension = Dimension(name='n')
+
         self._n_domains = kwargs.get('N', 1)
         self._global_bounds = kwargs.get('bounds', None)
-        self._implicit_exprs = None
 
     def __subdomain_finalize__(self, dimensions, shape, distributor=None, **kwargs):
         # Create the SubDomain's SubDimensions
-        sub_dimensions = []
-        for d in dimensions:
-            sub_dimensions.append(SubDimension.middle
-                                  ('%si_%s' % (d.name, self.implicit_dimension.name),
-                                   d, 0, 0))
-        self._dimensions = tuple(sub_dimensions)
+        self._dimensions = tuple(
+            MultiSubDimension('%si_%s' % (d.name, self.implicit_dimension.name), d)
+            for d in dimensions
+        )
 
         # Compute the SubDomainSet shapes
         global_bounds = []
@@ -621,57 +722,14 @@ class SubDomainSet(SubDomain):
         if distributor and distributor.is_parallel:
             # Now create local bounds based on distributor
             processed = []
-            for dim, d, m, M in zip(dimensions, distributor.decomposition, d_m, d_M):
-                bounds_m = np.zeros(m.shape, dtype=m.dtype)
-                bounds_M = np.zeros(m.shape, dtype=m.dtype)
-                for j in range(m.size):
-                    lmin = d.glb_min + m[j]
-                    lmax = d.glb_max - M[j]
-
-                    # Check if the subdomain doesn't intersect with the decomposition
-                    if lmin < d.loc_abs_min and lmax < d.loc_abs_min:
-                        bounds_m[j] = d.loc_abs_max
-                        bounds_M[j] = d.loc_abs_max
-                        continue
-                    if lmin > d.loc_abs_max and lmax > d.loc_abs_max:
-                        bounds_m[j] = d.loc_abs_max
-                        bounds_M[j] = d.loc_abs_max
-                        continue
-
-                    if lmin < d.loc_abs_min:
-                        bounds_m[j] = 0
-                    elif lmin > d.loc_abs_max:
-                        bounds_m[j] = d.loc_abs_max
-                        bounds_M[j] = d.loc_abs_max
-                        continue
-                    else:
-                        bounds_m[j] = d.index_glb_to_loc(m[j], LEFT)
-
-                    if lmax < d.loc_abs_min:
-                        bounds_m[j] = d.loc_abs_max
-                        bounds_M[j] = d.loc_abs_max
-                        continue
-                    elif lmax >= d.loc_abs_max:
-                        bounds_M[j] = 0
-                    else:
-                        bounds_M[j] = d.index_glb_to_loc(M[j], RIGHT)
-
-                processed.append(bounds_m)
-                processed.append(bounds_M)
+            for dec, m, M in zip(distributor.decomposition, d_m, d_M):
+                processed.extend(self._bounds_glb_to_loc(dec, m, M))
             self._local_bounds = as_tuple(processed)
         else:
             # Not distributed and hence local and global bounds are
             # equivalent.
             self._local_bounds = self._global_bounds
 
-    @property
-    def n_domains(self):
-        return self._n_domains
-
-    @property
-    def bounds(self):
-        return self._local_bounds
-
     def _create_implicit_exprs(self, grid):
         if not len(self._local_bounds) == 2*len(self.dimensions):
             raise ValueError("Left and right bounds must be supplied for each dimension")
@@ -697,3 +755,15 @@ class SubDomainSet(SubDomain):
                 func.data[:] = self._local_bounds[j]
             dat.append(Eq(d.thickness[j % 2][0], func[i_dim]))
         return as_tuple(dat)
+
+    @property
+    def n_domains(self):
+        return self._n_domains
+
+    @property
+    def bounds(self):
+        return self._local_bounds
+
+    @property
+    def implicit_dimensions(self):
+        return (self.implicit_dimension,)
diff --git a/tests/test_dse.py b/tests/test_dse.py
index 9ab538348..361f1a7b7 100644
--- a/tests/test_dse.py
+++ b/tests/test_dse.py
@@ -2565,7 +2565,7 @@ class TestTTI(object):
 
         # Check expected opcount/oi
         assert summary[('section1', None)].ops == 90
-        assert np.isclose(summary[('section1', None)].oi, 1.508, atol=0.001)
+        assert np.isclose(summary[('section1', None)].oi, 1.511, atol=0.001)
 
         # With optimizations enabled, there should be exactly four IncrDimensions
         op = wavesolver.op_fwd()
diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py
index 04a5e9ccf..dd83337ce 100644
--- a/tests/test_subdomains.py
+++ b/tests/test_subdomains.py
@@ -2,7 +2,9 @@ import pytest
 import numpy as np
 from math import floor
 
-from conftest import opts_tiling
+from sympy import sin, tan
+
+from conftest import opts_tiling, assert_structure
 from devito import (Grid, Function, TimeFunction, Eq, solve, Operator, SubDomain,
                     SubDomainSet, Dimension)
 from devito.ir import FindNodes, Expression
@@ -427,3 +429,69 @@ class TestSubdomains(object):
                              [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.int32)
 
         assert((np.array(f.data[:]+g.data[:]) == expected).all())
+
+    def test_issue_1761(self):
+        """
+        MFE for issue #1761.
+        """
+
+        class DummySubdomains(SubDomainSet):
+            name = 'dummydomain'
+        dummy = DummySubdomains(N=1, bounds=(1, 1, 1, 1))
+
+        grid = Grid(shape=(10, 10), subdomains=(dummy,))
+
+        f = TimeFunction(name='f', grid=grid)
+        g = TimeFunction(name='g', grid=grid)
+        theta = Function(name='theta', grid=grid)
+        phi = Function(name='phi', grid=grid)
+
+        eqns = [Eq(f.forward, f*sin(phi), subdomain=grid.subdomains['dummydomain']),
+                Eq(g.forward, g*sin(theta), subdomain=grid.subdomains['dummydomain'])]
+
+        op = Operator(eqns)
+
+        # Make sure it jit-compiles
+        op.cfunction
+
+        assert_structure(op, ['x,y', 't,n', 't,n,xi_n,yi_n'], 'x,y,t,n,xi_n,yi_n')
+
+    def test_issue_1761_b(self):
+        """
+        Follow-up issue emerged after patching #1761. The thicknesses assigments
+        were missing before the third equation.
+        """
+        n = Dimension(name='n')
+        m = Dimension(name='m')
+
+        class DummySubdomains(SubDomainSet):
+            name = 'dummydomain'
+            implicit_dimension = m
+
+        dummy = DummySubdomains(N=1, bounds=(1, 1, 1, 1))
+
+        class DummySubdomains2(SubDomainSet):
+            name = 'dummydomain2'
+            implicit_dimension = n
+
+        dummy2 = DummySubdomains2(N=1, bounds=(1, 1, 1, 1))
+
+        grid = Grid(shape=(10, 10), subdomains=(dummy, dummy2))
+
+        f = TimeFunction(name='f', grid=grid)
+        g = TimeFunction(name='g', grid=grid)
+        theta = Function(name='theta', grid=grid)
+        phi = Function(name='phi', grid=grid)
+
+        eqns = [Eq(f.forward, f*sin(phi), subdomain=grid.subdomains['dummydomain']),
+                Eq(g.forward, g*sin(theta), subdomain=grid.subdomains['dummydomain2']),
+                Eq(f.forward, f*tan(phi), subdomain=grid.subdomains['dummydomain'])]
+
+        op = Operator(eqns)
+
+        # Make sure it jit-compiles
+        op.cfunction
+
+        assert_structure(op, ['x,y', 't,m', 't,m,xi_m,yi_m',
+                              't,n', 't,n,xi_n,yi_n', 't,m', 't,m,xi_m,yi_m'],
+                         'x,y,t,m,xi_m,yi_m,n,xi_n,yi_n,m,xi_m,yi_m')
['devito/passes/clusters/misc.py', 'devito/ir/support/utils.py', 'devito/ir/clusters/algorithms.py', 'tests/test_subdomains.py', 'tests/test_dse.py', 'devito/core/cpu.py', 'devito/passes/clusters/aliases.py', 'devito/ir/support/space.py', 'devito/ir/equations/algorithms.py', 'devito/core/gpu.py', 'devito/types/grid.py']