/
domain_analysis.py
340 lines (276 loc) · 13.1 KB
/
domain_analysis.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
# -*- coding: utf-8 -*-
"""Algorithms for building canonical data structure for integrals over subdomains."""
# Copyright (C) 2009-2016 Anders Logg and Martin Sandve Alnæs
#
# This file is part of UFL (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
from collections import defaultdict
import ufl
from ufl.log import error
from ufl.integral import Integral
from ufl.form import Form
from ufl.sorting import cmp_expr, sorted_expr
from ufl.utils.sorting import canonicalize_metadata, sorted_by_key
from ufl.algorithms.coordinate_derivative_helpers import attach_coordinate_derivatives, strip_coordinate_derivatives
import numbers
class IntegralData(object):
"""Utility class with the members (domain, integral_type,
subdomain_id, integrals, metadata)
where metadata is an empty dictionary that may be used for
associating metadata with each object.
"""
__slots__ = ('domain', 'integral_type', 'subdomain_id',
'integrals', 'metadata',
'integral_coefficients',
'enabled_coefficients')
def __init__(self, domain, integral_type, subdomain_id, integrals,
metadata):
if 1 != len(set(itg.ufl_domain() for itg in integrals)):
error("Multiple domains mismatch in integral data.")
if not all(integral_type == itg.integral_type() for itg in integrals):
error("Integral type mismatch in integral data.")
if not all(subdomain_id == itg.subdomain_id() for itg in integrals):
error("Subdomain id mismatch in integral data.")
self.domain = domain
self.integral_type = integral_type
self.subdomain_id = subdomain_id
self.integrals = integrals
# This is populated in preprocess using data not available at
# this stage:
self.integral_coefficients = None
self.enabled_coefficients = None
# TODO: I think we can get rid of this with some refactoring
# in ffc:
self.metadata = metadata
def __lt__(self, other):
# To preserve behaviour of extract_integral_data:
return ((self.integral_type, self.subdomain_id,
self.integrals, self.metadata) <
(other.integral_type, other.subdomain_id, other.integrals,
other.metadata))
def __eq__(self, other):
# Currently only used for tests:
return (self.integral_type == other.integral_type and
self.subdomain_id == other.subdomain_id and
self.integrals == other.integrals and
self.metadata == other.metadata)
def __str__(self):
s = "IntegralData over domain(%s, %s), with integrals:\n%s\nand metadata:\n%s" % (
self.integral_type, self.subdomain_id,
'\n\n'.join(map(str, self.integrals)), self.metadata)
return s
def dicts_lt(a, b):
na = 0 if a is None else len(a)
nb = 0 if b is None else len(b)
if na != nb:
return len(a) < len(b)
for ia, ib in zip(sorted_by_key(a), sorted_by_key(b)):
# Assuming keys are sortable (usually str)
if ia[0] != ib[0]:
return (ia[0].__class__.__name__, ia[0]) < (ib[0].__class__.__name__, ib[0]) # Hack to preserve type sorting in py3
# Assuming values are sortable
if ia[1] != ib[1]:
return (ia[1].__class__.__name__, ia[1]) < (ib[1].__class__.__name__, ib[1]) # Hack to preserve type sorting in py3
# Tuple comparison helper
class ExprTupleKey(object):
__slots__ = ('x',)
def __init__(self, x):
self.x = x
def __lt__(self, other):
# Comparing expression first
c = cmp_expr(self.x[0], other.x[0])
if c < 0:
return True
elif c > 0:
return False
else:
# Comparing form compiler data
mds = canonicalize_metadata(self.x[1])
mdo = canonicalize_metadata(other.x[1])
return mds < mdo
def group_integrals_by_domain_and_type(integrals, domains):
"""
Input:
integrals: list of Integral objects
domains: list of AbstractDomain objects from the parent Form
Output:
integrals_by_domain_and_type: dict: (domain, integral_type) -> list(Integral)
"""
integrals_by_domain_and_type = defaultdict(list)
for itg in integrals:
if itg.ufl_domain() is None:
error("Integral has no domain.")
key = (itg.ufl_domain(), itg.integral_type())
# Append integral to list of integrals with shared key
integrals_by_domain_and_type[key].append(itg)
return integrals_by_domain_and_type
def integral_subdomain_ids(integral):
"Get a tuple of integer subdomains or a valid string subdomain from integral."
did = integral.subdomain_id()
if isinstance(did, numbers.Integral):
return (did,)
elif isinstance(did, tuple):
if not all(isinstance(d, numbers.Integral) for d in did):
error("Expecting only integer subdomains in tuple.")
return did
elif did in ("everywhere", "otherwise"):
# TODO: Define list of valid strings somewhere more central
return did
else:
error("Invalid domain id %s." % did)
def rearrange_integrals_by_single_subdomains(integrals, do_append_everywhere_integrals):
"""Rearrange integrals over multiple subdomains to single subdomain integrals.
Input:
integrals: list(Integral)
Output:
integrals: dict: subdomain_id -> list(Integral) (reconstructed with single subdomain_id)
"""
# Split integrals into lists of everywhere and subdomain integrals
everywhere_integrals = []
subdomain_integrals = []
for itg in integrals:
dids = integral_subdomain_ids(itg)
if dids == "otherwise":
error("'otherwise' integrals should never occur before preprocessing.")
elif dids == "everywhere":
everywhere_integrals.append(itg)
else:
subdomain_integrals.append((dids, itg))
# Fill single_subdomain_integrals with lists of integrals from
# subdomain_integrals, but split and restricted to single
# subdomain ids
single_subdomain_integrals = defaultdict(list)
for dids, itg in subdomain_integrals:
# Region or single subdomain id
for did in dids:
# Restrict integral to this subdomain!
single_subdomain_integrals[did].append(itg.reconstruct(subdomain_id=did))
# Add everywhere integrals to each single subdomain id integral
# list
otherwise_integrals = []
for ev_itg in everywhere_integrals:
# Restrict everywhere integral to 'otherwise'
otherwise_integrals.append(ev_itg.reconstruct(subdomain_id="otherwise"))
# Restrict everywhere integral to each subdomain
# and append to each integral list
if do_append_everywhere_integrals:
for subdomain_id in sorted(single_subdomain_integrals.keys()):
single_subdomain_integrals[subdomain_id].append(
ev_itg.reconstruct(subdomain_id=subdomain_id))
if otherwise_integrals:
single_subdomain_integrals["otherwise"] = otherwise_integrals
return single_subdomain_integrals
def accumulate_integrands_with_same_metadata(integrals):
"""
Taking input on the form:
integrals = [integral0, integral1, ...]
Return result on the form:
integrands_by_id = [(integrand0, metadata0),
(integrand1, metadata1), ...]
where integrand0 < integrand1 by the canonical ufl expression ordering criteria.
"""
# Group integrals by compiler data hash
by_cdid = {}
for itg in integrals:
cd = itg.metadata()
cdid = hash(canonicalize_metadata(cd))
if cdid not in by_cdid:
by_cdid[cdid] = ([], cd)
by_cdid[cdid][0].append(itg)
# Accumulate integrands separately for each compiler data object
# id
for cdid in by_cdid:
integrals, cd = by_cdid[cdid]
# Ensure canonical sorting of more than two integrands
integrands = sorted_expr((itg.integrand() for itg in integrals))
integrands_sum = sum(integrands[1:], integrands[0])
by_cdid[cdid] = (integrands_sum, cd)
# Sort integrands canonically by integrand first then compiler
# data
return sorted(by_cdid.values(), key=ExprTupleKey)
def build_integral_data(integrals):
"""Build integral data given a list of integrals.
:arg integrals: An iterable of :class:`~.Integral` objects.
:returns: A tuple of :class:`IntegralData` objects.
The integrals you pass in here must have been rearranged and
gathered (removing the "everywhere" subdomain_id. To do this, you
should call :func:`group_form_integrals`.
"""
itgs = defaultdict(list)
for integral in integrals:
domain = integral.ufl_domain()
integral_type = integral.integral_type()
subdomain_id = integral.subdomain_id()
if subdomain_id == "everywhere":
raise ValueError("'everywhere' not a valid subdomain id. Did you forget to call group_form_integrals?")
# Group for integral data (One integral data object for all
# integrals with same domain, itype, subdomain_id (but
# possibly different metadata).
itgs[(domain, integral_type, subdomain_id)].append(integral)
# Build list with canonical ordering, iteration over dicts
# is not deterministic across python versions
def keyfunc(item):
(d, itype, sid), integrals = item
return (d._ufl_sort_key_(), itype, (type(sid).__name__, sid))
integral_datas = []
for (d, itype, sid), integrals in sorted(itgs.items(), key=keyfunc):
integral_datas.append(IntegralData(d, itype, sid, integrals, {}))
return integral_datas
def group_form_integrals(form, domains, do_append_everywhere_integrals=True):
"""Group integrals by domain and type, performing canonical simplification.
:arg form: the :class:`~.Form` to group the integrals of.
:arg domains: an iterable of :class:`~.Domain`s.
:returns: A new :class:`~.Form` with gathered integrands.
"""
# Group integrals by domain and type
integrals_by_domain_and_type = \
group_integrals_by_domain_and_type(form.integrals(), domains)
integrals = []
for domain in domains:
for integral_type in ufl.measure.integral_types():
# Get integrals with this domain and type
ddt_integrals = integrals_by_domain_and_type.get((domain, integral_type))
if ddt_integrals is None:
continue
# Group integrals by subdomain id, after splitting e.g.
# f*dx((1,2)) + g*dx((2,3)) -> f*dx(1) + (f+g)*dx(2) + g*dx(3)
# (note: before this call, 'everywhere' is a valid subdomain_id,
# and after this call, 'otherwise' is a valid subdomain_id)
single_subdomain_integrals = \
rearrange_integrals_by_single_subdomains(ddt_integrals, do_append_everywhere_integrals)
for subdomain_id, ss_integrals in sorted_by_key(single_subdomain_integrals):
# strip the coordinate derivatives from all integrals
# this yields a list of the form [(coordinate derivative, integral), ...]
stripped_integrals_and_coordderivs = strip_coordinate_derivatives(ss_integrals)
# now group the integrals by the coordinate derivative
def calc_hash(cd):
return sum(sum(tuple_elem._ufl_compute_hash_()
for tuple_elem in tuple_) for tuple_ in cd)
coordderiv_integrals_dict = {}
for integral, coordderiv in stripped_integrals_and_coordderivs:
coordderivhash = calc_hash(coordderiv)
if coordderivhash in coordderiv_integrals_dict:
coordderiv_integrals_dict[coordderivhash][1].append(integral)
else:
coordderiv_integrals_dict[coordderivhash] = (coordderiv, [integral])
# cd_integrals_dict is now a dict of the form
# { hash: (CoordinateDerivative, [integral, integral, ...]), ... }
# we can now put the integrals back together and then afterwards
# apply the CoordinateDerivative again
for cdhash, samecd_integrals in sorted_by_key(coordderiv_integrals_dict):
# Accumulate integrands of integrals that share the
# same compiler data
integrands_and_cds = \
accumulate_integrands_with_same_metadata(samecd_integrals[1])
for integrand, metadata in integrands_and_cds:
integral = Integral(integrand, integral_type, domain,
subdomain_id, metadata, None)
integral = attach_coordinate_derivatives(integral, samecd_integrals[0])
integrals.append(integral)
return Form(integrals)
def reconstruct_form_from_integral_data(integral_data):
integrals = []
for ida in integral_data:
integrals.extend(ida.integrals)
return Form(integrals)