1
2
3 """Global function space discretization."""
4
5 from __future__ import division
6
7 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
8
9 __license__ = """
10 This program is free software: you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation, either version 3 of the License, or
13 (at your option) any later version.
14
15 This program is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with this program. If not, see U{http://www.gnu.org/licenses/}.
22 """
23
24
25
26
27 import numpy
28 import numpy.linalg as la
29 import hedge.tools
30 import hedge.mesh
31 import hedge.optemplate
32 import hedge._internal
33 from pytools import memoize_method
34
35
36
37
38 -class _FaceGroup(hedge._internal.FaceGroup):
39 - def __init__(self, double_sided, debug):
43
45 return self.fil_registry.register(identifier, generator)
46
47 - def commit(self, discr, ldis_loc, ldis_opp):
85
90 """Once fully filled, this structure has the following data members:
91
92 @ivar members: a list of L{hedge.mesh.Element} instances in this group.
93 @ivar member_nrs: a list of the element ID numbers in this group.
94 @ivar local_discretization: an instance of L{hedge.element.Element}.
95 @ivar ranges: a list of C{slice} objects indicating the DOF numbers for
96 each element. Note: This is actually a C++ ElementRanges object.
97 @ivar mass_matrix: The element-local mass matrix M{M}.
98 @ivar inverse_mass_matrix: the element-local inverese mass matrix M{M^{-1}}.
99 @ivar differentiation_matrices: local differentiation matrices M{D_r, D_s, D_t},
100 i.e. differentiation by M{r, s, t, ....}.
101 @ivar stiffness_matrices: the element-local stiffness matrices M{M*D_r, M*D_s,...}.
102 @ivar jacobians: list of jacobians over all elements
103 @ivar inverse_jacobians: inverses of L{jacobians}.
104 @ivar diff_coefficients: a M{(d,d)}-matrix of coefficient vectors to turn
105 M{(r,s,t)}-differentiation into M{(x,y,z)}.
106 """
107 pass
108
113 - def __init__(self, nodes, ranges, vol_indices, face_groups,
114 el_face_to_face_group_and_face_pair={}):
115 self.nodes = nodes
116 self.ranges = ranges
117 self.vol_indices = vol_indices
118 self.face_groups = face_groups
119 self.el_face_to_face_group_and_face_pair = \
120 el_face_to_face_group_and_face_pair
121
123 fg, fp_idx = self.el_face_to_face_group_and_face_pair[el_face]
124
125 return fg.face_pairs[fp_idx]
126
128 fp = self.find_facepair(el_face)
129 el, face_nbr = el_face
130
131 for flux_face in [fp.loc, fp.opp]:
132 if flux_face.element_id == el.id and flux_face.face_id == face_nbr:
133 return flux_face
134 raise KeyError, "flux face not found in boundary"
135
140 - def __init__(self, discr, pp_optemplate):
141 self.discr = discr
142 self.pp_optemplate = pp_optemplate
143
145 return self.discr.run_preprocessed_optemplate(self.pp_optemplate, vars)
146
151 - def __init__(self, discr, el_range, interp_coeff):
152 self.discr = discr
153 self.el_range = el_range
154 self.interp_coeff = interp_coeff
155
157 from hedge.tools import log_shape
158 ls = log_shape(field)
159 if ls != ():
160 result = numpy.zeros(ls, dtype=self.discr.default_scalar_type)
161 from pytools import indices_in_shape
162 for i in indices_in_shape(ls):
163 result[i] = numpy.dot(
164 self.interp_coeff, field[i][self.el_range])
165 return result
166 else:
167 return numpy.dot(self.interp_coeff, field[self.el_range])
168
171 """The global approximation space.
172
173 Instances of this class tie together a local discretization (i.e. polynomials
174 on an elemnent) into a function space on a mesh. They provide creation
175 functions such as interpolating given functions, differential operators and
176 flux lifting operators.
177 """
178
179 @classmethod
181 return set([
182 "ilist_generation",
183 "node_permutation",
184 "print_op_code",
185 "dump_dataflow_graph",
186 ])
187
188 @classmethod
190 return set([
191 "ilist_generation",
192 "node_permutation",
193 ])
194
195 @staticmethod
197 if local_discretization is None and order is None:
198 raise ValueError, "must supply either local_discretization or order"
199 if local_discretization is not None and order is not None:
200 raise ValueError, "must supply only one of local_discretization and order"
201 if local_discretization is None:
202 from hedge.element import ELEMENTS
203 from pytools import one
204 ldis_class = one(
205 ldis_class for ldis_class in ELEMENTS
206 if isinstance(mesh.elements[0], ldis_class.geometry))
207 return ldis_class(order)
208 else:
209 return local_discretization
210
211 - def __init__(self, mesh, local_discretization=None,
212 order=None, debug=set(), default_scalar_type=numpy.float64,
213 run_context=None):
241
244
245
247 from pytools.log import IntervalTimer
248
249 self.gather_timer = IntervalTimer("t_gather",
250 "Time spent gathering fluxes")
251 self.lift_timer = IntervalTimer("t_lift",
252 "Time spent lifting fluxes")
253 self.mass_timer = IntervalTimer("t_mass",
254 "Time spent applying mass operators")
255 self.diff_timer = IntervalTimer("t_diff",
256 "Time spent applying applying differentiation operators")
257 self.vector_math_timer = IntervalTimer("t_vector_math",
258 "Time spent doing vector math")
259
260 return [self.gather_timer,
261 self.lift_timer,
262 self.mass_timer,
263 self.diff_timer,
264 self.vector_math_timer]
265
267 from pytools.log import IntervalTimer, EventCounter
268
269 self.gather_counter = EventCounter("n_gather",
270 "Number of flux gather invocations")
271 self.lift_counter = EventCounter("n_lift",
272 "Number of flux lift invocations")
273 self.mass_counter = EventCounter("n_mass_op",
274 "Number of mass operator applications")
275 self.diff_counter = EventCounter("n_diff",
276 "Number of differentiation operator applications")
277
278 self.gather_flop_counter = EventCounter("n_flops_gather",
279 "Number of floating point operations in gather")
280 self.lift_flop_counter = EventCounter("n_flops_lift",
281 "Number of floating point operations in lift")
282 self.mass_flop_counter = EventCounter("n_flops_mass",
283 "Number of floating point operations in mass operator")
284 self.diff_flop_counter = EventCounter("n_flops_diff",
285 "Number of floating point operations in diff operator")
286 self.vector_math_flop_counter = EventCounter("n_flops_vector_math",
287 "Number of floating point operations in vector math")
288
289 self.interpolant_counter = EventCounter("n_interp",
290 "Number of interpolant evaluations")
291
292 self.interpolant_timer = IntervalTimer("t_interp",
293 "Time spent evaluating interpolants")
294
295 for op in self.create_op_timers():
296 mgr.add_quantity(op)
297
298 mgr.add_quantity(self.gather_counter)
299 mgr.add_quantity(self.lift_counter)
300 mgr.add_quantity(self.mass_counter)
301 mgr.add_quantity(self.diff_counter)
302
303 mgr.add_quantity(self.gather_flop_counter)
304 mgr.add_quantity(self.lift_flop_counter)
305 mgr.add_quantity(self.mass_flop_counter)
306 mgr.add_quantity(self.diff_flop_counter)
307 mgr.add_quantity(self.vector_math_flop_counter)
308
309 mgr.add_quantity(self.interpolant_counter)
310 mgr.add_quantity(self.interpolant_timer)
311
312 from pytools.log import time_and_count_function
313 self.interpolate_volume_function = \
314 time_and_count_function(
315 self.interpolate_volume_function,
316 self.interpolant_timer,
317 self.interpolant_counter)
318
319 self.interpolate_boundary_function = \
320 time_and_count_function(
321 self.interpolate_boundary_function,
322 self.interpolant_timer,
323 self.interpolant_counter)
324
325 from pytools import single_valued
326 try:
327 order = single_valued(eg.local_discretization.order
328 for eg in self.element_groups)
329 except ValueError:
330 pass
331 else:
332 mgr.set_constant("dg_order", order)
333
334 mgr.set_constant("default_type", self.default_scalar_type.__name__)
335 mgr.set_constant("element_count", len(self.mesh.elements))
336 mgr.set_constant("node_count", len(self.nodes))
337
338 for f in self.all_debug_flags():
339 mgr.set_constant("debug_%s" % f, f in self.debug)
340
341 self.instrumented = True
342
343
345 from hedge._internal import UniformElementRanges
346
347 eg = _ElementGroup()
348 eg.members = self.mesh.elements
349 eg.member_nrs = numpy.fromiter((el.id for el in eg.members), dtype=numpy.uint32)
350 eg.local_discretization = ldis = local_discretization
351 eg.ranges = UniformElementRanges(
352 0,
353 len(ldis.unit_nodes()),
354 len(self.mesh.elements))
355
356 nodes_per_el = ldis.node_count()
357
358
359
360
361
362
363
364
365
366 self.nodes = numpy.empty(
367 (len(self.mesh.elements)*nodes_per_el, self.dimensions),
368 dtype=float, order="C")
369
370 unit_nodes = numpy.empty( (nodes_per_el, self.dimensions),
371 dtype=float, order="C")
372
373 for i_node, node in enumerate(ldis.unit_nodes()):
374 unit_nodes[i_node] = node
375
376 from hedge._internal import map_element_nodes
377
378 for el in self.mesh.elements:
379 map_element_nodes(
380 self.nodes,
381 el.id*nodes_per_el*self.dimensions,
382 el.map,
383 unit_nodes,
384 self.dimensions)
385
386 self.group_map = [(eg, i) for i in range(len(self.mesh.elements))]
387 self.element_groups = [eg]
388
390 for eg in self.element_groups:
391 ldis = eg.local_discretization
392
393 mmat = eg.mass_matrix = ldis.mass_matrix()
394 immat = eg.inverse_mass_matrix = ldis.inverse_mass_matrix()
395 dmats = eg.differentiation_matrices = \
396 ldis.differentiation_matrices()
397 smats = eg.stiffness_matrices = \
398 [numpy.dot(mmat, d) for d in dmats]
399 smats = eg.stiffness_t_matrices = \
400 [numpy.dot(d.T, mmat.T) for d in dmats]
401 eg.minv_st = \
402 [numpy.dot(numpy.dot(immat,d.T), mmat) for d in dmats]
403
404 eg.jacobians = numpy.array([
405 abs(el.map.jacobian())
406 for el in eg.members])
407 eg.inverse_jacobians = numpy.array([
408 abs(el.inverse_map.jacobian())
409 for el in eg.members])
410
411 eg.diff_coefficients = numpy.array([
412 [
413 [
414 el.inverse_map
415 .matrix[loc_coord, glob_coord]
416 for el in eg.members
417 ]
418 for loc_coord in range(ldis.dimensions)
419 ]
420 for glob_coord in range(ldis.dimensions)
421 ])
422
423 eg.stiffness_coefficients = numpy.array([
424 [
425 [
426 abs(el.map.jacobian())*el.inverse_map
427 .matrix[loc_coord, glob_coord]
428 for el in eg.members
429 ]
430 for loc_coord in range(ldis.dimensions)
431 ]
432 for glob_coord in range(ldis.dimensions)
433 ])
434
449
451 from hedge._internal import FacePair
452 from hedge.element import FaceVertexMismatch
453
454 fg = _FaceGroup(double_sided=True,
455 debug="ilist_generation" in self.debug)
456
457 all_ldis_l = []
458 all_ldis_n = []
459
460 debug_node_perm = "node_permutation" in self.debug
461
462
463 for i, (local_face, neigh_face) in enumerate(self.mesh.interfaces):
464 e_l, fi_l = local_face
465 e_n, fi_n = neigh_face
466
467 eslice_l, ldis_l = self.find_el_data(e_l.id)
468 eslice_n, ldis_n = self.find_el_data(e_n.id)
469
470 all_ldis_l.append(ldis_l)
471 all_ldis_n.append(ldis_n)
472
473 vertices_l = e_l.faces[fi_l]
474 vertices_n = e_n.faces[fi_n]
475
476 findices_l = ldis_l.face_indices()[fi_l]
477 findices_n = ldis_n.face_indices()[fi_n]
478
479 try:
480 findices_shuffle_op_n = \
481 ldis_l.get_face_index_shuffle_to_match(
482 vertices_l, vertices_n)
483
484 if debug_node_perm and ldis_l.has_facial_nodes and ldis_n.has_facial_nodes:
485 findices_shuffled_n = findices_shuffle_op_n(findices_n)
486
487 for i, j in zip(findices_l, findices_shuffled_n):
488 dist = self.nodes[eslice_l.start+i]-self.nodes[eslice_n.start+j]
489 assert la.norm(dist) < 1e-14
490
491 except FaceVertexMismatch:
492
493
494
495 vertices_n, axis = self.mesh.periodic_opposite_faces[vertices_n]
496
497 findices_shuffle_op_n = \
498 ldis_l.get_face_index_shuffle_to_match(vertices_l, vertices_n)
499
500 if debug_node_perm and ldis_l.has_facial_nodes and ldis_n.has_facial_nodes:
501 findices_shuffled_n = findices_shuffle_op_n(findices_n)
502
503 for i, j in zip(findices_l, findices_shuffled_n):
504 dist = self.nodes[eslice_l.start+i]-self.nodes[eslice_n.start+j]
505 dist[axis] = 0
506 assert la.norm(dist) < 1e-14
507
508
509 fp = FacePair()
510
511 fp.loc.el_base_index = eslice_l.start
512 fp.opp.el_base_index = eslice_n.start
513
514 fp.loc.face_index_list_number = fg.register_face_index_list(
515 identifier=fi_l,
516 generator=lambda: findices_l)
517 fp.opp.face_index_list_number = fg.register_face_index_list(
518 identifier=(fi_n, findices_shuffle_op_n),
519 generator=lambda : findices_shuffle_op_n(findices_n))
520 from pytools import get_write_to_map_from_permutation
521 fp.opp_native_write_map = fg.register_face_index_list(
522 identifier=(fi_n, findices_shuffle_op_n, "wtm"),
523 generator=lambda :
524 get_write_to_map_from_permutation(
525 findices_shuffle_op_n(findices_n), findices_n))
526
527 self._set_flux_face_data(fp.loc, ldis_l, local_face)
528 self._set_flux_face_data(fp.opp, ldis_n, neigh_face)
529
530
531 fp.loc.h = fp.opp.h = max(fp.loc.h, fp.opp.h)
532
533 assert len(fp.__dict__) == 0
534 assert len(fp.loc.__dict__) == 0
535 assert len(fp.opp.__dict__) == 0
536
537 fg.face_pairs.append(fp)
538
539 if len(fg.face_pairs):
540 from pytools import single_valued
541 ldis_l = single_valued(all_ldis_l)
542 ldis_n = single_valued(all_ldis_n)
543
544 fg.commit(self, ldis_l, ldis_n)
545
546 self.face_groups = [fg]
547 else:
548 self.face_groups = []
549
551 return bool(self.mesh.tag_to_boundary.get(tag, []))
552
553 @memoize_method
555 """Get a _Boundary instance for a given `tag'.
556
557 If there is no boundary tagged with `tag', an empty _Boundary instance
558 is returned. Asking for a nonexistant boundary is not an error.
559 (Otherwise get_boundary would unnecessarily become non-local when run
560 in parallel.)
561 """
562 from hedge._internal import FacePair
563
564 nodes = []
565 face_ranges = {}
566 vol_indices = []
567 face_group = _FaceGroup(double_sided=False,
568 debug="ilist_generation" in self.debug)
569 ldis = None
570 el_face_to_face_group_and_face_pair = {}
571
572 for ef in self.mesh.tag_to_boundary.get(tag, []):
573 el, face_nr = ef
574
575 el_slice, ldis = self.find_el_data(el.id)
576 face_indices = ldis.face_indices()[face_nr]
577
578 f_start = len(nodes)
579 nodes += [self.nodes[el_slice.start+i] for i in face_indices]
580 face_ranges[ef] = (f_start, len(nodes))
581 vol_indices.extend(el_slice.start+i for i in face_indices)
582
583
584 fp = FacePair()
585 fp.loc.el_base_index = el_slice.start
586 fp.opp.el_base_index = f_start
587 fp.loc.face_index_list_number = face_group.register_face_index_list(
588 identifier=face_nr,
589 generator=lambda: face_indices)
590 fp.opp.face_index_list_number = face_group.register_face_index_list(
591 identifier=(),
592 generator=lambda: tuple(xrange(len(face_indices))))
593 self._set_flux_face_data(fp.loc, ldis, ef)
594 assert len(fp.__dict__) == 0
595 assert len(fp.loc.__dict__) == 0
596 assert len(fp.opp.__dict__) == 0
597
598 face_group.face_pairs.append(fp)
599
600
601 el_face_to_face_group_and_face_pair[ef] = \
602 face_group, len(face_group.face_pairs)-1
603
604 face_group.commit(self, ldis, ldis)
605
606 bdry = _Boundary(
607 nodes=numpy.array(nodes),
608 ranges=face_ranges,
609 vol_indices=numpy.asarray(vol_indices, dtype=numpy.intp),
610 face_groups=[face_group],
611 el_face_to_face_group_and_face_pair=
612 el_face_to_face_group_and_face_pair)
613
614 return bdry
615
616
618 """Return the number of nodes in this discretization."""
619 return len(self.nodes)
620
623
626
627 compute_kind = "numpy"
628
630 from hedge.tools import with_object_array_or_scalar
631 if dtype is not None:
632 return with_object_array_or_scalar(lambda f: f.astype(dtype), field)
633 else:
634 return field
635
637 orig_kind = self.get_kind(field)
638
639 if orig_kind != "numpy":
640 raise ValueError, "unable to perform kind conversion: %s -> %s" % (
641 orig_kind, kind)
642
643 return field
644
646 orig_kind = self.get_kind(field)
647
648 if orig_kind != "numpy":
649 raise ValueError, "unable to perform kind conversion: %s -> %s" % (
650 orig_kind, kind)
651
652 return field
653
671
672 - def volume_empty(self, shape=(), dtype=None, kind="numpy"):
673 if kind != "numpy":
674 raise ValueError, "invalid vector kind requested"
675
676 if dtype is None:
677 dtype = self.default_scalar_type
678 return numpy.empty(shape+(len(self.nodes),), dtype)
679
680 - def volume_zeros(self, shape=(), dtype=None, kind="numpy"):
681 if kind != "numpy":
682 raise ValueError, "invalid vector kind requested"
683
684 if dtype is None:
685 dtype = self.default_scalar_type
686 return numpy.zeros(shape+(len(self.nodes),), dtype)
687
689 if kind is None:
690 kind = self.compute_kind
691
692 try:
693
694 shape = f.shape
695 except AttributeError:
696
697 shape = ()
698
699 slice_pfx = (slice(None),)*len(shape)
700 out = self.volume_empty(shape, dtype, kind="numpy")
701 for eg in self.element_groups:
702 for el, el_slice in zip(eg.members, eg.ranges):
703 for point_nr in xrange(el_slice.start, el_slice.stop):
704 out[slice_pfx + (point_nr,)] = \
705 f(self.nodes[point_nr], el)
706 return self.convert_volume(out, kind=kind)
707
709 if kind not in ["numpy", "numpy-mpi-recv"]:
710 raise ValueError, "invalid vector kind requested"
711
712 if dtype is None:
713 dtype = self.default_scalar_type
714 return numpy.empty(shape+(len(self.get_boundary(tag).nodes),), dtype)
715
717 if kind not in ["numpy", "numpy-mpi-recv"]:
718 raise ValueError, "invalid vector kind requested"
719 if dtype is None:
720 dtype = self.default_scalar_type
721
722 return numpy.zeros(shape+(len(self.get_boundary(tag).nodes),), dtype)
723
741
742 @memoize_method
757
759 if kind is None:
760 kind = self.compute_kind
761
762 if kind != "numpy":
763 raise ValueError("invalid target vector kind in volumize_boundary_field")
764
765 bdry = self.get_boundary(tag)
766
767 def f(subfld):
768 result = self.volume_zeros(dtype=bfield.dtype, kind="numpy")
769 result[bdry.vol_indices] = subfld
770 return result
771
772 from hedge.tools import with_object_array_or_scalar
773 return with_object_array_or_scalar(f, bfield)
774
776 if kind is None:
777 kind = self.compute_kind
778
779 if kind != "numpy":
780 raise ValueError("invalid target vector kind in boundarize_volume_field")
781
782 bdry = self.get_boundary(tag)
783
784 from hedge.tools import log_shape, is_obj_array
785 ls = log_shape(field)
786
787 if is_obj_array(field):
788 if len(field) == 0:
789 return numpy.zeros(())
790
791 result = self.boundary_empty(tag, shape=ls, dtype=field[0].dtype)
792 from pytools import indices_in_shape
793 for i in indices_in_shape(ls):
794 result[i] = field[i][bdry.vol_indices]
795
796 return result
797 else:
798 return field[tuple(slice(None) for i in range(len(ls))) + (bdry.vol_indices,)]
799
804
806 return numpy.array(indices, dtype=numpy.intp)
807
808
810 return numpy.dot(a, b)
811
812 @memoize_method
816
817 @memoize_method
820
822 from hedge.tools import log_shape
823
824 ls = log_shape(volume_vector)
825 if ls == ():
826 if isinstance(volume_vector, (int, float, complex)):
827
828 empty = self.volume_empty(dtype=type(volume_vector))
829 empty.fill(volume_vector)
830 volume_vector = empty
831
832 return self.nodewise_dot_product(
833 self._mass_ones(), volume_vector)
834 else:
835 result = numpy.zeros(shape=ls, dtype=float)
836
837 from pytools import indices_in_shape
838 for i in indices_in_shape(ls):
839 vvi = volume_vector[i]
840 if isinstance(vvi, (int, float)) and vvi == 0:
841 result[i] = 0
842 else:
843 result[i] = self.nodewise_dot_product(
844 self._mass_ones(), volume_vector[i])
845
846 return result
847
848 @memoize_method
853
854 - def norm(self, volume_vector, p=2):
855 if p == numpy.Inf:
856 return numpy.abs(volume_vector).max()
857 else:
858 from hedge.tools import log_shape
859
860 if p != 2:
861 volume_vector = numpy.abs(volume_vector)**(p/2)
862
863 return self.inner_product(
864 volume_vector,
865 volume_vector)**(1/p)
866
868 mass_op = self._compiled_mass_operator()
869
870 from hedge.tools import log_shape
871 ls = log_shape(a)
872 assert log_shape(b) == ls
873 if ls == ():
874 return float(self.nodewise_dot_product(
875 a, mass_op(b)))
876 else:
877 assert len(ls) == 1
878 return float(sum(
879 self.nodewise_dot_product(
880 sub_a, mass_op(sub_b))
881 for sub_a, sub_b in zip(a,b)))
882
883
885 group, idx = self.group_map[el_id]
886 return group.ranges[idx]
887
889 return self.group_map[el_id][0].local_discretization
890
892 group, idx = self.group_map[el_id]
893 return group.ranges[idx], group.local_discretization
894
896 for i, (start, stop) in enumerate(self.element_group):
897 if start <= idx < stop:
898 return i
899 raise ValueError, "not a valid dof index"
900
901
903 distinct_ldis = set(eg.local_discretization for eg in self.element_groups)
904 return min(ldis.dt_non_geometric_factor()
905 for ldis in distinct_ldis)
906
912
913 - def dt_factor(self, max_system_ev, stepper_class=None, *stepper_args):
980
982 for eg in self.element_groups:
983 for el, rng in zip(eg.members, eg.ranges):
984 if el.contains_point(point):
985 ldis = eg.local_discretization
986 basis_values = numpy.array([
987 phi(el.inverse_map(point))
988 for phi in ldis.basis_functions()])
989 vdm_t = ldis.vandermonde().T
990 return _PointEvaluator(
991 discr=self,
992 el_range=rng,
993 interp_coeff=la.solve(vdm_t, basis_values))
994
995 raise RuntimeError, "point %s not found" % point
996
997
998 - def compile(self, optemplate, post_bind_mapper=lambda x: x):
999 ex = self.executor_class(self, optemplate, post_bind_mapper)
1000
1001 if "dump_dataflow_graph" in self.debug:
1002 ex.code.dump_dataflow_graph()
1003
1004 if self.instrumented:
1005 ex.instrument()
1006 return ex
1007
1009 self.exec_functions[name] = func
1010
1016 """A symmetry map on global DG functions.
1017
1018 Suppose that the L{Mesh} on which a L{Discretization} is defined has
1019 is mapped onto itself by a nontrivial symmetry map M{f(.)}. Then
1020 this class allows you to carry out this map on vectors representing
1021 functions on this L{Discretization}.
1022 """
1023 - def __init__(self, discr, sym_map, element_map, threshold=1e-13):
1024 self.discretization = discr
1025
1026 complete_el_map = {}
1027 for i, j in element_map.iteritems():
1028 complete_el_map[i] = j
1029 complete_el_map[j] = i
1030
1031 self.map = {}
1032
1033 for eg in discr.element_groups:
1034 for el, el_slice in zip(eg.members, eg.ranges):
1035 mapped_i_el = complete_el_map[el.id]
1036 mapped_slice = discr.find_el_range(mapped_i_el)
1037 for i_pt in range(el_slice.start, el_slice.stop):
1038 pt = discr.nodes[i_pt]
1039 mapped_pt = sym_map(pt)
1040 for m_i_pt in range(mapped_slice.start, mapped_slice.stop):
1041 if la.norm(discr.nodes[m_i_pt] - mapped_pt) < threshold:
1042 self.map[i_pt] = m_i_pt
1043 break
1044
1045 if i_pt not in self.map:
1046 for m_i_pt in range(mapped_slice.start, mapped_slice.stop):
1047 print la.norm_2(discr.nodes[m_i_pt] - mapped_pt)
1048 raise RuntimeError, "no symmetry match found"
1049
1055
1060 result = discr.volume_zeros()
1061 import random
1062 for eg in discr.element_groups:
1063 for e_start, e_end in eg.ranges:
1064 result[e_start:e_end] = random.random()
1065 return result
1066
1071 result = discr.volume_zeros()
1072
1073 try:
1074 faces = discr.mesh.tag_to_boundary[tag]
1075 except KeyError:
1076 pass
1077 else:
1078 for face in faces:
1079 el, fl = face
1080
1081 (el_start, el_end), ldis = discr.find_el_data(el.id)
1082 fl_indices = ldis.face_indices()[fl]
1083
1084 for i in fl_indices:
1085 result[el_start+i] = 1
1086
1087 return result
1088
1093 result = discr.volume_empty()
1094 result.fill(1)
1095 return result
1096
1102 - def __init__(self, from_discr, to_discr):
1103 self.from_discr = from_discr
1104 self.to_discr = to_discr
1105
1106 self.interp_matrices = []
1107 for from_eg, to_eg in zip(
1108 from_discr.element_groups, to_discr.element_groups):
1109 from_ldis = from_eg.local_discretization
1110 to_ldis = to_eg.local_discretization
1111
1112 from_count = from_ldis.node_count()
1113 to_count = to_ldis.node_count()
1114
1115
1116 for from_el, to_el in zip(from_eg.members, to_eg.members):
1117 assert from_el is to_el
1118
1119 from hedge.tools import permutation_matrix
1120
1121
1122
1123 if to_count > from_count:
1124 to_node_ids_to_idx = dict(
1125 (nid, i) for i, nid in
1126 enumerate(to_ldis.generate_mode_identifiers()))
1127
1128 to_indices = [
1129 to_node_ids_to_idx[from_nid]
1130 for from_nid in from_ldis.generate_mode_identifiers()
1131 ]
1132
1133 pmat = permutation_matrix(
1134 to_indices=to_indices,
1135 h=to_count, w=from_count)
1136 else:
1137 from_node_ids_to_idx = dict(
1138 (nid, i) for i, nid in
1139 enumerate(from_ldis.generate_mode_identifiers()))
1140
1141 from_indices = [
1142 from_node_ids_to_idx[to_nid]
1143 for to_nid in to_ldis.generate_mode_identifiers()
1144 ]
1145
1146 pmat = permutation_matrix(
1147 from_indices=from_indices,
1148 h=to_count, w=from_count)
1149
1150
1151 from_matrix = from_ldis.vandermonde()
1152 to_matrix = to_ldis.vandermonde()
1153
1154 from hedge.tools import leftsolve
1155 from numpy import dot
1156 self.interp_matrices.append(
1157 numpy.asarray(
1158 leftsolve(from_matrix, dot(to_matrix, pmat)),
1159 order="C"))
1160
1162 from hedge._internal import perform_elwise_operator
1163 from hedge.tools import log_shape
1164
1165 ls = log_shape(from_vec)
1166 result = self.to_discr.volume_zeros(ls, kind="numpy")
1167
1168 from pytools import indices_in_shape
1169 for i in indices_in_shape(ls):
1170 for from_eg, to_eg, imat in zip(
1171 self.from_discr.element_groups,
1172 self.to_discr.element_groups,
1173 self.interp_matrices):
1174 perform_elwise_operator(
1175 from_eg.ranges, to_eg.ranges,
1176 imat, from_vec[i], result[i])
1177
1178 return result
1179
1185 """A typical exponential-falloff mode response filter function.
1186
1187 See description in Section 5.6.1 of Hesthaven/Warburton.
1188 """
1189 - def __init__(self, min_amplification=0.1, order=6):
1190 """Construct the filter function.
1191
1192 The amplification factor of the lowest-order (constant) mode is always 1.
1193
1194 @arg min_amplification: The amplification factor applied to the highest mode.
1195 @arg order: The order of the filter. This controls how fast (or slowly) the
1196 C{min_amplification} is reached.
1197 """
1198 from math import log
1199 self.alpha = -log(min_amplification)
1200 self.order = order
1201
1203 eta = sum(mode_idx)/ldis.order
1204
1205 from math import exp
1206 return exp(-self.alpha * eta**self.order)
1207
1212 - def __init__(self, discr, mode_response_func):
1213 """Construct a filter.
1214
1215 @arg discr: The L{Discretization} for which the filter is to be
1216 constructed.
1217 @arg mode_response_func: A function mapping
1218 C{(mode_tuple, local_discretization)} to a float indicating the
1219 factor by which this mode is to be multiplied after filtering.
1220 """
1221 self.discr = discr
1222
1223 self.filter_map = {}
1224
1225 for eg in discr.element_groups:
1226 ldis = eg.local_discretization
1227
1228 node_count = ldis.node_count()
1229
1230 filter_coeffs = [mode_response_func(mid, ldis)
1231 for mid in ldis.generate_mode_identifiers()]
1232
1233
1234 vdm = ldis.vandermonde()
1235 from hedge.tools import leftsolve
1236 from numpy import dot
1237 mat = numpy.asarray(
1238 leftsolve(vdm,
1239 dot(vdm, numpy.diag(filter_coeffs))),
1240 order="C")
1241 self.filter_map[eg] = mat
1242
1244 from hedge.tools import log_shape
1245
1246 ls = log_shape(vec)
1247 result = self.discr.volume_zeros(ls)
1248
1249 from pytools import indices_in_shape
1250 for i in indices_in_shape(ls):
1251 from hedge._internal import perform_elwise_operator
1252 for eg in self.discr.element_groups:
1253 perform_elwise_operator(eg.ranges, eg.ranges,
1254 self.filter_map[eg], vec[i], result[i])
1255
1256 return result
1257
1259 return self.filter_map[el_group]
1260