gistfile1.diff

=== modified file 'python/fluidity/bcs.py'
--- python/fluidity/bcs.py	2013-05-29 13:25:20 +0000
+++ python/fluidity/bcs.py	2013-05-31 09:14:46 +0000
@@ -82,39 +82,54 @@
         # Add method argument if it's given
         if "method" in kwargs:
             args = tuple(list(args) + [kwargs["method"]])
+        self._cached_setup = {}
+
+    def _setup(self, domain, b):
+        rank = 0 if b.dim == (1,) else len(b.dim)
+        key = (domain, rank, b.cdim)
+        try:
+            # elem_map, domain_nodes, kernel
+            return self._cached_setup[key]
+        except KeyError:
+            pass
+        mesh = self.field.mesh
+        # Get list of nodes in domain.
+        # FIXME: get zero_rows to take a map
+        # so we don't mandate pulling data back and forth in
+        # the interface.
+        domain_nodes = set()
+        # If running in parallel, this processor's mesh partition might
+        # not contain any nodes of the current domain
+        if domain not in mesh.faces.boundary_list:
+            assert op2.MPI.parallel, \
+                "Domain %d not in the boundary list for mesh %s" % (domain, mesh)
+            elem_map = op2.Map(op2.Set(0), mesh.node_set(rank),
+                               mesh.faces.surface_mesh.shape.loc,
+                               numpy.zeros(0))
+        else:
+            elem_map = mesh.faces.surface_elements_to_nodes_maps[domain][rank]
+            for nodelist in elem_map.values:
+                for node in nodelist:
+                    if node < mesh.node_classes[1]:
+                        for i in range(b.cdim):
+                            domain_nodes.add(b.cdim * node + i)
+        domain_nodes = list(domain_nodes)
+        kernel = self._apply_dirichlet(elem_map.dim, b.cdim)
+        self._cached_setup[key] = elem_map, domain_nodes, kernel
+        return self._cached_setup[key]
 
     def apply(self, A, b):
         rank = 0 if b.dim == (1,) else len(b.dim)
         mesh = self.field.mesh
         mesh.compute_boundaries()
         for domain in self.domains:
-            # Get list of nodes in domain.
-            # FIXME: get zero_rows to take a map
-            # so we don't mandate pulling data back and forth in
-            # the interface.
-            domain_nodes = set()
-            # If running in parallel, this processor's mesh partition might
-            # not contain any nodes of the current domain
-            if domain not in mesh.faces.boundary_list:
-                assert op2.MPI.parallel, \
-                    "Domain %d not in the boundary list for mesh %s" % (domain, mesh)
-                elem_map = op2.Map(op2.Set(0), mesh.node_set(rank),
-                                   mesh.faces.surface_mesh.shape.loc,
-                                   numpy.zeros(0))
-            else:
-                elem_map = mesh.faces.surface_elements_to_nodes_maps[domain][rank]
-                for nodelist in elem_map.values:
-                    for node in nodelist:
-                        if node < mesh.node_classes[1]:
-                            for i in range(b.cdim):
-                                domain_nodes.add(b.cdim * node + i)
-            domain_nodes = list(domain_nodes)
+            elem_map, domain_nodes, kernel = self._setup(domain, b)
             A.zero_rows(domain_nodes, 1.0)
             # Apply BC on RHS
             # FIXME: Need to generate the appropriate code for setting the BC depending
             # on the number of nodes and data dimension
             value = op2.Global(b.dim, self.expression)
-            op2.par_loop(self._apply_dirichlet(elem_map.dim, b.cdim),
+            op2.par_loop(kernel,
                          elem_map.iterset,
                          b(elem_map, op2.WRITE),
                          value(op2.READ))