[analyticalproblems] use LincombFunction to represent affine decompos…

…itions

src/pymor/analyticalproblems/elliptic.py

@@ -11,13 +11,11 @@
 
 
 class EllipticProblem(ImmutableInterface):
-    """Affinely decomposed linear elliptic problem.
+    """Linear elliptic problem description.
 
     The problem consists in solving ::
 
-    |        Kd                                     Kv                             Kr
-    | - ∇ ⋅  ∑  θ_{d,k}(μ) ⋅ d_k(x) ∇ u(x, μ) + ∇ ⋅ ∑  θ_{v,k}(μ) v_k(x) u(x, μ) + ∑  θ_{r,k}(μ) r_k(x) u(x, μ) = f(x, μ)
-    |       k=0                                    k=0                            k=0
+        - ∇ ⋅ [d(x, μ) ∇ u(x, μ)] + ∇ ⋅ [v(x, μ) u(x, μ)] + c(x, μ) u(x, μ) = f(x, μ)
 
     for u.
 
@@ -28,27 +26,13 @@ class EllipticProblem(ImmutableInterface):
     rhs
         The |Function| f(x, μ). `rhs.dim_domain` has to agree with the
         dimension of `domain`, whereas `rhs.shape_range` has to be `()`.
-    diffusion_functions
-        List containing the |Functions| d_k(x), each having `shape_range`
-        of either `()` or `(dim domain, dim domain)`.
-    diffusion_functionals
-        List containing the |ParameterFunctionals| θ_{d,k}(μ). If
-        `len(diffusion_functions) == 1`, `diffusion_functionals` is allowed
-        to be `None`, in which case no parameter dependence is assumed.
-    advection_functions
-        List containing the |Functions| v_k(x), each having `shape_range`
-        of `(dim domain,)`.
-    advection_functionals
-        List containing the |ParameterFunctionals| θ_{v,k}(μ). If
-        `len(advection_functions) == 1`, `advection_functionals` is allowed
-        to be `None`, in which case no parameter dependence is assumed.
-    reaction_functions
-        List containing the |Functions| r_k(x), each having `shape_range`
-        of `()`.
-    reaction_functionals
-        List containing the |ParameterFunctionals| θ_{r,k}(μ). If
-        `len(reaction_functions) == 1`, `reaction_functionals` is allowed
-        to be `None`, in which case no parameter dependence is assumed.
+    diffusion
+        The |Function| d(x, μ) with `shape_range` of either `()` or
+        `(dim domain, dim domain)`.
+    advection
+        The |Function| v(x, μ), with `shape_range` of `(dim domain,)`.
+    reaction
+        The |Function| c(x, μ), with `shape_range` of `()`.
     dirichlet_data
         |Function| providing the Dirichlet boundary values.
     neumann_data
@@ -64,71 +48,31 @@ class EllipticProblem(ImmutableInterface):
     ----------
     domain
     rhs
-    diffusion_functions
-    diffusion_functionals
-    advection_functions
-    advection_functionals
-    reaction_functions
-    reaction_functionals
+    diffusion
+    advection
+    reaction
     dirichlet_data
     neumann_data
     robin_data
     """
 
     def __init__(self, domain=RectDomain(), rhs=ConstantFunction(dim_domain=2),
-                 diffusion_functions=None,
-                 diffusion_functionals=None,
-                 advection_functions=None,
-                 advection_functionals=None,
-                 reaction_functions=None,
-                 reaction_functionals=None,
+                 diffusion=None, advection=None, reaction=None,
                  dirichlet_data=None, neumann_data=None, robin_data=None,
                  parameter_space=None, name=None):
-        assert diffusion_functions is None or isinstance(diffusion_functions, (tuple, list))
-        assert advection_functions is None or isinstance(advection_functions, (tuple, list))
-        assert reaction_functions is None or isinstance(reaction_functions, (tuple, list))
-
-        assert diffusion_functionals is None and diffusion_functions is None \
-            or diffusion_functionals is None and len(diffusion_functions) == 1 \
-            or len(diffusion_functionals) == len(diffusion_functions)
-        assert advection_functionals is None and advection_functions is None \
-            or advection_functionals is None and len(advection_functions) == 1 \
-            or len(advection_functionals) == len(advection_functions)
-        assert reaction_functionals is None and reaction_functions is None \
-            or reaction_functionals is None and len(reaction_functions) == 1 \
-            or len(reaction_functionals) == len(reaction_functions)
-
-        # for backward compatibility:
-        if (diffusion_functions is None and advection_functions is None and reaction_functions is None):
-            diffusion_functions = (ConstantFunction(dim_domain=2),)
-
-        # dim_domain:
-        if diffusion_functions is not None:
-            dim_domain = diffusion_functions[0].dim_domain
-
-        assert rhs.dim_domain == dim_domain
-        if diffusion_functions is not None:
-            for f in diffusion_functions:
-                assert f.dim_domain == dim_domain
-        if advection_functions is not None:
-            for f in advection_functions:
-                assert f.dim_domain == dim_domain
-        if reaction_functions is not None:
-            for f in reaction_functions:
-                assert f.dim_domain == dim_domain
-
-        assert dirichlet_data is None or dirichlet_data.dim_domain == dim_domain
-        assert neumann_data is None or neumann_data.dim_domain == dim_domain
+        assert rhs.dim_domain == domain.dim
+        assert diffusion is None or diffusion.dim_domain == domain.dim
+        assert advection is None or advection.dim_domain == domain.dim
+        assert reaction is None or reaction.dim_domain == domain.dim
+        assert dirichlet_data is None or dirichlet_data.dim_domain == domain.dim
+        assert neumann_data is None or neumann_data.dim_domain == domain.dim
         assert robin_data is None or (isinstance(robin_data, tuple) and len(robin_data) == 2)
-        assert robin_data is None or np.all([f.dim_domain == dim_domain for f in robin_data])
+        assert robin_data is None or np.all([f.dim_domain == domain.dim for f in robin_data])
         self.domain = domain
         self.rhs = rhs
-        self.diffusion_functions = diffusion_functions
-        self.diffusion_functionals = diffusion_functionals
-        self.advection_functions = advection_functions
-        self.advection_functionals = advection_functionals
-        self.reaction_functions = reaction_functions
-        self.reaction_functionals = reaction_functionals
+        self.diffusion = diffusion
+        self.advection = advection
+        self.reaction = reaction
         self.dirichlet_data = dirichlet_data
         self.neumann_data = neumann_data
         self.robin_data = robin_data

src/pymor/analyticalproblems/helmholtz.py

@@ -5,7 +5,7 @@
 
 from pymor.analyticalproblems.elliptic import EllipticProblem
 from pymor.domaindescriptions.basic import RectDomain
-from pymor.functions.basic import ConstantFunction
+from pymor.functions.basic import ConstantFunction, LincombFunction
 from pymor.parameters.functionals import ExpressionParameterFunctional
 from pymor.parameters.spaces import CubicParameterSpace
 
@@ -46,13 +46,10 @@ def helmholtz_problem(domain=RectDomain(), rhs=None, parameter_range=(0., 100.),
 
         neumann_data=neumann_data,
 
-        diffusion_functions=[ConstantFunction(1., dim_domain=domain.dim)],
+        diffusion=ConstantFunction(1., dim_domain=domain.dim),
 
-        diffusion_functionals=[1.],
-
-        reaction_functions=[ConstantFunction(1., dim_domain=domain.dim)],
-
-        reaction_functionals=[ExpressionParameterFunctional('-k**2', {'k': ()})],
+        reaction=LincombFunction([ConstantFunction(1., dim_domain=domain.dim)],
+                                 [ExpressionParameterFunctional('-k**2', {'k': ()})]),
 
         parameter_space=CubicParameterSpace({'k': ()}, *parameter_range),
 

src/pymor/analyticalproblems/thermalblock.py

@@ -7,7 +7,7 @@
 
 from pymor.analyticalproblems.elliptic import EllipticProblem
 from pymor.domaindescriptions.basic import RectDomain
-from pymor.functions.basic import ConstantFunction, ExpressionFunction
+from pymor.functions.basic import ConstantFunction, ExpressionFunction, LincombFunction
 from pymor.parameters.functionals import ProjectionParameterFunctional
 from pymor.parameters.spaces import CubicParameterSpace
 
@@ -67,11 +67,11 @@ def diffusion_function_factory(ix, iy):
 
         rhs=ConstantFunction(dim_domain=2, value=1.),
 
-        diffusion_functions=[diffusion_function_factory(ix, iy)
-                             for ix, iy in product(range(num_blocks[0]), range(num_blocks[1]))],
-
-        diffusion_functionals=[parameter_functional_factory(ix, iy)
-                               for ix, iy in product(range(num_blocks[0]), range(num_blocks[1]))],
+        diffusion=LincombFunction([diffusion_function_factory(ix, iy)
+                                   for ix, iy in product(range(num_blocks[0]), range(num_blocks[1]))],
+                                  [parameter_functional_factory(ix, iy)
+                                   for ix, iy in product(range(num_blocks[0]), range(num_blocks[1]))],
+                                  name='diffusion'),
 
         parameter_space=CubicParameterSpace({'diffusion': (num_blocks[1], num_blocks[0])}, *parameter_range),
 

src/pymor/discretizers/elliptic.py

@@ -5,6 +5,7 @@
 from pymor.analyticalproblems.elliptic import EllipticProblem
 from pymor.discretizations.basic import StationaryDiscretization
 from pymor.domaindiscretizers.default import discretize_domain_default
+from pymor.functions.basic import LincombFunction
 from pymor.grids.boundaryinfos import EmptyBoundaryInfo
 from pymor.grids.referenceelements import line, triangle, square
 from pymor.gui.qt import PatchVisualizer, Matplotlib1DVisualizer
@@ -71,88 +72,48 @@ def discretize_elliptic_cg(analytical_problem, diameter=None, domain_discretizer
 
     p = analytical_problem
 
-    if p.diffusion_functionals is not None or p.advection_functionals is not None or p.reaction_functionals is not None:
-        # parametric case
-        Li = [DiffusionOperator(grid, boundary_info, diffusion_constant=0, name='boundary_part')]
-        coefficients = [1.]
-
-        # diffusion part
-        if p.diffusion_functionals is not None:
-            Li += [DiffusionOperator(grid, boundary_info, diffusion_function=df, dirichlet_clear_diag=True,
-                                     name='diffusion_{}'.format(i))
-                   for i, df in enumerate(p.diffusion_functions)]
-            coefficients += list(p.diffusion_functionals)
-        elif p.diffusion_functions is not None:
-            assert len(p.diffusion_functions) == 1
-            Li += [DiffusionOperator(grid, boundary_info, diffusion_function=p.diffusion_functions[0],
-                                     dirichlet_clear_diag=True, name='diffusion')]
-            coefficients.append(1.)
-
-        # advection part
-        if p.advection_functionals is not None:
-            Li += [AdvectionOperator(grid, boundary_info, advection_function=af, dirichlet_clear_diag=True,
-                                     name='advection_{}'.format(i))
-                   for i, af in enumerate(p.advection_functions)]
-            coefficients += list(p.advection_functionals)
-        elif p.advection_functions is not None:
-            assert len(p.advection_functions) == 1
-            Li += [AdvectionOperator(grid, boundary_info, advection_function=p.advection_functions[0],
-                                     dirichlet_clear_diag=True, name='advection')]
-            coefficients.append(1.)
-
-        # reaction part
-        if p.reaction_functionals is not None:
-            Li += [ReactionOperator(grid, boundary_info, coefficient_function=rf, dirichlet_clear_diag=True,
-                                    name='reaction_{}'.format(i))
-                   for i, rf in enumerate(p.reaction_functions)]
-            coefficients += list(p.reaction_functionals)
-        elif p.reaction_functions is not None:
-            assert len(p.reaction_functions) == 1
-            Li += [ReactionOperator(grid, boundary_info, coefficient_function=p.reaction_functions[0],
-                                    dirichlet_clear_diag=True, name='reaction')]
-            coefficients.append(1.)
-
-        # robin boundaries
-        if p.robin_data is not None:
-            Li += [cg.RobinBoundaryOperator(grid, boundary_info, robin_data=p.robin_data, order=2, name='robin')]
-            coefficients.append(1.)
-
-        L = LincombOperator(operators=Li, coefficients=coefficients, name='ellipticOperator')
-    else:
-        # unparametric case, not operator for boundary treatment
-        Li = []
-
-        # only one operator has diagonal values, all subsequent operators have clear_diag
-        dirichlet_clear_diag = False
-        # diffusion part
-        if p.diffusion_functions is not None:
-            assert len(p.diffusion_functions) == 1
-            Li += [DiffusionOperator(grid, boundary_info, diffusion_function=p.diffusion_functions[0],
-                                     dirichlet_clear_diag=dirichlet_clear_diag, name='diffusion')]
-            dirichlet_clear_diag = True
-
-        # advection part
-        if p.advection_functions is not None:
-            assert len(p.advection_functions) == 1
-            Li += [AdvectionOperator(grid, boundary_info, advection_function=p.advection_functions[0],
-                                     dirichlet_clear_diag=dirichlet_clear_diag, name='advection')]
-            dirichlet_clear_diag = True
-
-        # reaction part
-        if p.reaction_functions is not None:
-            assert len(p.reaction_functions) == 1
-            Li += [ReactionOperator(grid, boundary_info, coefficient_function=p.reaction_functions[0],
-                                    dirichlet_clear_diag=dirichlet_clear_diag, name='reaction')]
-            dirichlet_clear_diag = True
-
-        # robin boundaries
-        if p.robin_data is not None:
-            Li += [cg.RobinBoundaryOperator(grid, boundary_info, robin_data=p.robin_data, order=2, name='robin')]
-
-        if len(Li) == 1:
-            L = Li[0]
-        else:
-            L = LincombOperator(operators=Li, coefficients=[1.] * len(Li), name='ellipticOperator')
+    Li = [DiffusionOperator(grid, boundary_info, diffusion_constant=0, name='boundary_part')]
+    coefficients = [1.]
+
+    # diffusion part
+    if isinstance(p.diffusion, LincombFunction):
+        Li += [DiffusionOperator(grid, boundary_info, diffusion_function=df, dirichlet_clear_diag=True,
+                                 name='diffusion_{}'.format(i))
+               for i, df in enumerate(p.diffusion.functions)]
+        coefficients += list(p.diffusion.coefficients)
+    elif p.diffusion is not None:
+        Li += [DiffusionOperator(grid, boundary_info, diffusion_function=p.diffusion,
+                                 dirichlet_clear_diag=True, name='diffusion')]
+        coefficients.append(1.)
+
+    # advection part
+    if isinstance(p.advection, LincombFunction):
+        Li += [AdvectionOperator(grid, boundary_info, advection_function=af, dirichlet_clear_diag=True,
+                                 name='advection_{}'.format(i))
+               for i, af in enumerate(p.advection.functions)]
+        coefficients += list(p.advection.coefficients)
+    elif p.advection is not None:
+        Li += [AdvectionOperator(grid, boundary_info, advection_function=p.advection,
+                                 dirichlet_clear_diag=True, name='advection')]
+        coefficients.append(1.)
+
+    # reaction part
+    if isinstance(p.reaction, LincombFunction):
+        Li += [ReactionOperator(grid, boundary_info, coefficient_function=rf, dirichlet_clear_diag=True,
+                                name='reaction_{}'.format(i))
+               for i, rf in enumerate(p.reaction.functions)]
+        coefficients += list(p.reaction.coefficients)
+    elif p.reaction is not None:
+        Li += [ReactionOperator(grid, boundary_info, coefficient_function=p.reaction,
+                                dirichlet_clear_diag=True, name='reaction')]
+        coefficients.append(1.)
+
+    # robin boundaries
+    if p.robin_data is not None:
+        Li += [cg.RobinBoundaryOperator(grid, boundary_info, robin_data=p.robin_data, order=2, name='robin')]
+        coefficients.append(1.)
+
+    L = LincombOperator(operators=Li, coefficients=coefficients, name='ellipticOperator')
 
     F = Functional(grid, p.rhs, boundary_info, dirichlet_data=p.dirichlet_data, neumann_data=p.neumann_data)
 
@@ -221,9 +182,9 @@ def discretize_elliptic_fv(analytical_problem, diameter=None, domain_discretizer
     assert boundary_info is None or grid is not None
     assert grid is None or domain_discretizer is None
 
-    if analytical_problem.advection_functions is not None:
+    if analytical_problem.advection is not None:
         raise NotImplementedError
-    if analytical_problem.reaction_functions is not None:
+    if analytical_problem.reaction is not None:
         raise NotImplementedError
     if analytical_problem.robin_data is not None:
         raise NotImplementedError
@@ -237,29 +198,27 @@ def discretize_elliptic_fv(analytical_problem, diameter=None, domain_discretizer
 
     p = analytical_problem
 
-    if p.diffusion_functionals is not None:
+    if isinstance(p.diffusion, LincombFunction):
         Li = [fv.DiffusionOperator(grid, boundary_info, diffusion_function=df, name='diffusion_{}'.format(i))
-              for i, df in enumerate(p.diffusion_functions)]
-        L = LincombOperator(operators=Li, coefficients=list(p.diffusion_functionals),
+              for i, df in enumerate(p.diffusion.functions)]
+        L = LincombOperator(operators=Li, coefficients=list(p.diffusion.coefficients),
                             name='diffusion')
 
         F0 = fv.L2ProductFunctional(grid, p.rhs, boundary_info=boundary_info, neumann_data=p.neumann_data)
         if p.dirichlet_data is not None:
             Fi = [fv.L2ProductFunctional(grid, None, boundary_info=boundary_info, dirichlet_data=p.dirichlet_data,
                                          diffusion_function=df, name='dirichlet_{}'.format(i))
-                  for i, df in enumerate(p.diffusion_functions)]
-            F = LincombOperator(operators=[F0] + Fi, coefficients=[1.] + list(p.diffusion_functionals),
+                  for i, df in enumerate(p.diffusion.functions)]
+            F = LincombOperator(operators=[F0] + Fi, coefficients=[1.] + list(p.diffusion.coefficients),
                                 name='rhs')
         else:
             F = F0
 
     else:
-        assert len(p.diffusion_functions) == 1
-        L = fv.DiffusionOperator(grid, boundary_info, diffusion_function=p.diffusion_functions[0],
-                                 name='diffusion')
+        L = fv.DiffusionOperator(grid, boundary_info, diffusion_function=p.diffusion, name='diffusion')
 
         F = fv.L2ProductFunctional(grid, p.rhs, boundary_info=boundary_info, dirichlet_data=p.dirichlet_data,
-                                   diffusion_function=p.diffusion_functions[0], neumann_data=p.neumann_data)
+                                   diffusion_function=p.diffusion, neumann_data=p.neumann_data)
 
     if grid.reference_element in (triangle, square):
         visualizer = PatchVisualizer(grid=grid, bounding_box=grid.bounding_box(), codim=0)