Commits

Martin Alnæs committed 9b9329c

Work on quadrature rule code generation.

  • Participants
  • Parent commits cc7e914

Comments (0)

Files changed (1)

sandbox/integralgeneratordesign/integralgenerator.py

 
+from uflacs.utils.log import debug, info, warning, error, uflacs_assert
 from uflacs.codeutils.format_code_structure import format_code_structure, Indented, Block, ForRange
 from uflacs.geometry.default_names import names
 
+# TODO: Refactor
+from uflacs.backends.ffc.ffc_statement_formatter import CppStatementFormatterRules
+langfmt = CppStatementFormatterRules()
+
 
 class IntegralGenerator(object):
     def __init__(self):
         # TODO: This is just fake test data
-
-        self._argument_space_dimensions = [7, 8] # ir["prim_idims"] # FIXME: *2 for dS?
-        rank = len(self._argument_space_dimensions)
-        self._idofs = ["%s%d" % (names.ia, i) for i in range(rank)]
-
-
-        self._num_points = [1, 2]
-        self._dof_ranges = {
+        ir = {}
+        ir["prim_idims"] = [7, 8]
+        ir["quadrature_weights"] = {
+            1: ( [0.5, 0.5],  [(0.1, 0.2), (0.3, 0.4)] ),
+            2: ( [3.5, 3.5],  [(1.1, 1.2), (2.3, 2.4)] ),
+            }
+        ir["dof_ranges"] = {
             1: [
                 ((0,3), (3,6)),
                 ((1,2), (4,5))
                 ((0,5),)
                 ],
             }
-        self._num_arguments = 2
+
+        # ...
+
+        quadrature_rules = ir["quadrature_weights"]
+        self._num_points = []
+        self._weights = {}
+        self._points = {}
+        for num_points in sorted(quadrature_rules.keys()):
+            self._num_points.append(num_points)
+            w, p = quadrature_rules[num_points]
+            self._weights[num_points] = w
+            self._points[num_points] = p
+
+        if len(quadrature_rules) > 1:
+            warning("Multiple quadrature rules not fully implemented, will likely crash somewhere.")
+
+        self._argument_space_dimensions = ir["prim_idims"] # FIXME: *2 for dS?
+        self._dof_ranges = ir["dof_ranges"]
+
+        self._num_arguments = len(self._argument_space_dimensions)
+        self._idofs = ["%s%d" % (names.ia, i) for i in range(self._num_arguments)]
+
 
     def generate(self):
         """Generate entire tabulate_tensor body."""
         parts = []
+        parts += [self.generate_static_tables()]
         parts += [self.generate_pre_quadrature()]
         parts += [self.generate_quadrature_loops()]
         parts += [self.generate_post_quadrature()]
         return format_code_structure(Indented(parts))
 
+    def generate_static_tables(self): # FIXME
+        """Generate static tables of precomputed values.
+
+        This typically includes quadrature points and weights,
+        as well as basis function values in quadrature points.
+        """
+        parts = ["// Quadrature weights and points"]
+
+        for num_points in self._num_points:
+            weights = self._weights[num_points]
+            points = self._points[num_points]
+            pdim = len(points[0])
+
+            weights = "{ %s }" % langfmt.precision_floats(weights)
+            points = "{ %s }" % langfmt.precision_floats(x for p in points for x in p)
+
+            wname = "%s%d" % (names.weights, num_points) # FIXME: Make sure we use this everywhere
+            pname = "%s%d" % (names.points, num_points)
+
+            parts += [langfmt.array_decl("static const double", wname, num_points, weights)]
+            parts += [langfmt.array_decl("static const double", pname, num_points*pdim, points)]
+            parts += [""]
+
+        return parts
+
     def generate_pre_quadrature(self): # FIXME
         """Generate initialization statements.