/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py

"""This module provides functionality for form assembly in Python,
corresponding to the C++ assembly and PDE classes.

The C++ :py:class:`assemble <dolfin.cpp.assemble>` function
(renamed to cpp_assemble) is wrapped with an additional
preprocessing step where code is generated using the
FFC JIT compiler.

The C++ PDE classes are reimplemented in Python since the C++ classes
rely on the dolfin::Form class which is not used on the Python side."""

# Copyright (C) 2007-2008 Anders Logg
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
#
# Modified by Martin Sandve Alnaes, 2008.
# Modified by Johan Hake, 2008-2009.
# Modified by Garth N. Wells, 2008-2013.
# Modified by Joachim B. Haga, 2012.
#
# First added:  2007-08-15
# Last changed: 2013-02-18

__all__ = ["assemble", "assemble_system", "SystemAssembler"]

import types

# UFL modules
from ufl import interval, triangle, tetrahedron

# Import SWIG-generated extension module (DOLFIN C++)
import dolfin.cpp as cpp

# Local imports
from dolfin.fem.form import *

# JIT assembler
def assemble(form,
             tensor=None,
             mesh=None,
             coefficients=None,
             function_spaces=None,
             cell_domains=None,
             exterior_facet_domains=None,
             interior_facet_domains=None,
             reset_sparsity=True,
             add_values=False,
             finalize_tensor=True,
             keep_diagonal=False,
             backend=None,
             form_compiler_parameters=None,
             bcs=None):
    """
    Assemble the given form and return the corresponding tensor.

    *Arguments*
        Depending on the input form, which may be a functional, linear
        form, bilinear form or higher rank form, a scalar value, a vector,
        a matrix or a higher rank tensor is returned.

    In the simplest case, no additional arguments are needed. However,
    additional arguments may and must in some cases be provided as
    outlined below.

    The ``form`` can be either an FFC form or a precompiled UFC
    form. If a precompiled or 'pure' UFC form is given, then
    ``coefficients`` and ``function_spaces`` have to be provided
    too. The coefficient functions should be provided as a 'dict'
    using the FFC functions as keys. The function spaces should be
    provided either as a list where the number of function spaces must
    correspond to the number of basis functions in the form, or as a
    single argument, implying that the same FunctionSpace is used for
    all test/trial spaces.

    If the form defines integrals over different subdomains,
    :py:class:`MeshFunctions <dolfin.cpp.MeshFunction>` over the
    corresponding topological entities defining the subdomains can be
    provided. An instance of a :py:class:`SubDomain
    <dolfin.cpp.SubDomain>` can also be passed for each subdomain.

    The implementation of the returned tensor is determined by the
    default linear algebra backend. This can be overridden by
    specifying a different backend.

    Each call to assemble() will create a new tensor. If the
    ``tensor`` argument is provided, this will be used instead. If
    ``reset_sparsity`` is set to False, the provided tensor will not be
    reset to zero before assembling (adding) values to the tensor.

    If the ``keep_diagonal`` is set to True, assembler ensures that
    potential zeros on a matrix diagonal are kept in sparsity pattern
    so every diagonal entry can be changed in a future (for example
    by ident() or ident_zeros()).

    Specific form compiler parameters can be provided by the
    ``form_compiler_parameters`` argument. Form compiler parameters
    can also be controlled using the global parameters stored in
    parameters["form_compiler"].

    *Examples of usage*
        The standard stiffness matrix ``A`` and a load vector ``b``
        can be assembled as follows:

        .. code-block:: python

            A = assemble(inner(grad(u),grad(v))*dx())
            b = assemble(f*v*dx())

        It is possible to explicitly prescribe the domains over which
        integrals wll be evaluated using the arguments
        ``cell_domains``, ``exterior_facet_domains`` and
        ``interior_facet_domains``. For instance, using a mesh
        function marking parts of the boundary:

        .. code-block:: python

            # MeshFunction marking boundary parts
            boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

            # Sample variational forms
            a = inner(grad(u), grad(v))*dx() + p*u*v*ds(0)
            L = f*v*dx() - g*v*ds(1) + p*q*v*ds(0)

            A = assemble(a, exterior_facet_domains=boundary_parts)
            b = assemble(L, exterior_facet_domains=boundary_parts)

        To ensure that the assembled matrix has the right type, one may use
        the ``tensor`` argument:

        .. code-block:: python

            A = PETScMatrix()
            assemble(a, tensor=A)

        The form ``a`` is now assembled into the PETScMatrix ``A``.

    """

    # Extract common cell from mesh (may be missing in form definition)
    common_cell = None if mesh is None else mesh.ufl_cell()

    # First check if we got a cpp.Form which originates from cpp layer
    if isinstance(form, cpp.Form) and not hasattr(form, "_compiled_form"):

        # Then we just try to use that one
        dolfin_form = form

    else:

        # Wrap form
        dolfin_form = Form(form,
                           function_spaces=function_spaces,
                           coefficients=coefficients,
                           form_compiler_parameters=form_compiler_parameters,
                           common_cell=common_cell)

    # Set mesh if specified (important for functionals without a function spaces)
    if mesh is not None:
        dolfin_form.set_mesh(mesh)

    # Create tensor
    tensor = _create_tensor(form, dolfin_form.rank(), backend, tensor)

    # Extract domains
    cell_domains, exterior_facet_domains, interior_facet_domains = \
        _extract_domains(dolfin_form.mesh(),
                         cell_domains,
                         exterior_facet_domains,
                         interior_facet_domains)

    # Attach domains to form 
    if cell_domains is not None:
        dolfin_form.set_cell_domains(cell_domains);
    if interior_facet_domains is not None:
        dolfin_form.set_interior_facet_domains(interior_facet_domains);
    if exterior_facet_domains is not None:
        dolfin_form.set_exterior_facet_domains(exterior_facet_domains);

    # Call C++ assemble function
    assembler = cpp.Assembler()
    assembler.reset_sparsity = reset_sparsity
    assembler.add_values = add_values
    assembler.finalize_tensor = finalize_tensor
    assembler.keep_diagonal = keep_diagonal
    assembler.assemble(tensor,
                 dolfin_form)

    # Convert to float for scalars
    if dolfin_form.rank() == 0:
        tensor = tensor.getval()

    # Apply (possibly list of) boundary conditions
    bcs = _wrap_in_list(bcs, 'bcs', cpp.DirichletBC)
    if bcs:
        cpp.deprecation("Passing \"bcs\" to assemble", "1.3.0",
                        "Apply DirichletBC manually after an assemble.")
    for bc in bcs:
        bc.apply(tensor)

    # Return value
    return tensor

# JIT system assembler
def assemble_system(A_form,
                    b_form,
                    bcs=None,
                    x0=None,
                    A_coefficients=None,
                    b_coefficients=None,
                    A_function_spaces=None,
                    b_function_spaces=None,
                    cell_domains=None,
                    exterior_facet_domains=None,
                    interior_facet_domains=None,
                    reset_sparsity=True,
                    add_values=False,
                    finalize_tensor=True,
                    keep_diagonal=False,
                    A_tensor=None,
                    b_tensor=None,
                    backend=None,
                    form_compiler_parameters=None):
    """
    Assemble form(s) and apply any given boundary conditions in a
    symmetric fashion and return tensor(s).

    The standard application of boundary conditions does not
    necessarily preserve the symmetry of the assembled matrix. In
    order to perserve symmetry in a system of equations with boundary
    conditions, one may use the alternative assemble_system instead of
    multiple calls to :py:func:`assemble
    <dolfin.fem.assembling.assemble>`.

    *Examples of usage*

       For instance, the statements

       .. code-block:: python

           A = assemble(a)
           b = assemble(L)
           bc.apply(A, b)

       can alternatively be carried out by

       .. code-block:: python

           A, b = assemble_system(a, L, bc)

       The statement above is valid even if ``bc`` is a list of
       :py:class:`DirichletBC <dolfin.fem.bcs.DirichletBC>`
       instances. For more info and options, see :py:func:`assemble
       <dolfin.fem.assembling.assemble>`.

    """

    # Extract subdomains
    subdomains = { "cell": cell_domains,
                   "exterior_facet": exterior_facet_domains,
                   "interior_facet": interior_facet_domains}

    # Wrap forms
    A_dolfin_form = Form(A_form, A_function_spaces, A_coefficients,
                         subdomains, form_compiler_parameters)
    b_dolfin_form = Form(b_form, b_function_spaces, b_coefficients,
                         subdomains, form_compiler_parameters)

    # Create tensors
    A_tensor = _create_tensor(A_form, A_dolfin_form.rank(), backend, A_tensor)
    b_tensor = _create_tensor(b_form, b_dolfin_form.rank(), backend, b_tensor)

    # Extract domains
    cell_domains, exterior_facet_domains, interior_facet_domains = \
        _extract_domains(A_dolfin_form.mesh(),
                         cell_domains,
                         exterior_facet_domains,
                         interior_facet_domains)

    # Attach domains to form (from Mesh, in other cases form already
    # has domain data attached)
    if cell_domains is not None:
        A_dolfin_form.set_cell_domains(cell_domains);
        b_dolfin_form.set_cell_domains(cell_domains);
    if interior_facet_domains is not None:
        A_dolfin_form.set_interior_facet_domains(interior_facet_domains);
        b_dolfin_form.set_interior_facet_domains(interior_facet_domains);
    if exterior_facet_domains is not None:
        A_dolfin_form.set_exterior_facet_domains(exterior_facet_domains);
        b_dolfin_form.set_exterior_facet_domains(exterior_facet_domains);

    # Check bcs
    bcs = _wrap_in_list(bcs, 'bcs', cpp.DirichletBC)

    # Call C++ assemble function
    assembler = cpp.SystemAssembler(A_dolfin_form, b_dolfin_form, bcs)
    assembler.reset_sparsity = reset_sparsity
    assembler.add_values = add_values
    assembler.finalize_tensor = finalize_tensor
    assembler.keep_diagonal = keep_diagonal
    if x0 is not None:
        assembler.assemble(A_tensor, b_tensor, x0)
    else:
        assembler.assemble(A_tensor, b_tensor)

    return A_tensor, b_tensor

def _wrap_in_list(obj, name, types=type):
    if obj is None:
        lst = []
    elif hasattr(obj, '__iter__'):
        lst = list(obj)
    else:
        lst = [obj]
    for obj in lst:
        if not isinstance(obj, types):
            raise TypeError("expected a (list of) %s as '%s' argument" % (str(types),name))
    return lst

def _create_tensor(form, rank, backend, tensor):
    "Create tensor for form"

    # Check if tensor is supplied by user
    if tensor is not None:
        return tensor

    # Check backend argument
    if (not backend is None) and (not isinstance(backend, cpp.GenericLinearAlgebraFactory)):
        raise TypeError, "Provide a GenericLinearAlgebraFactory as 'backend'"

    # Create tensor
    if rank == 0:
        tensor = cpp.Scalar()
    elif rank == 1:
        if backend: tensor = backend.create_vector()
        else:       tensor = cpp.Vector()
    elif rank == 2:
        if backend: tensor = backend.create_matrix()
        else:       tensor = cpp.Matrix()
    else:
        raise RuntimeError, "Unable to create tensors of rank %d." % rank

    return tensor

def _extract_domains(mesh,
                     cell_domains,
                     exterior_facet_domains,
                     interior_facet_domains):

    def check_domain_type(domain, domain_type):
        if not isinstance(domain, (cpp.SubDomain,cpp.MeshFunctionSizet, types.NoneType)):
            raise TypeError, "expected a 'SubDomain', 'MeshFunction' of 'Sizet' or 'None', for the '%s'"%domain_type

    def build_mf(subdomain, dim):
        " Builds a MeshFunction from a SubDomain"
        mf = cpp.MeshFunction("size_t", mesh, dim)
        mf.set_all(1)
        subdomain.mark(mf,0)
        return mf

    # Type check of input
    check_domain_type(cell_domains, "cell_domains")
    check_domain_type(exterior_facet_domains, "exterior_facet_domains")
    check_domain_type(interior_facet_domains, "interior_facet_domains")

    # The cell dimension
    cell_dim = mesh.topology().dim()

    # Get cell_domains
    if isinstance(cell_domains, cpp.SubDomain):
        cell_domains = build_mf(cell_domains, cell_dim)

    # Get exterior_facet_domains (may be stored as part of the mesh)
    if isinstance(exterior_facet_domains, cpp.SubDomain):
        exterior_facet_domains = build_mf(exterior_facet_domains, cell_dim-1)

    # Get interior_facet_domains
    if isinstance(interior_facet_domains, cpp.SubDomain):
        interior_facet_domains = build_mf(interior_facet_domains, cell_dim-1)

    return cell_domains, exterior_facet_domains, interior_facet_domains

class SystemAssembler(cpp.SystemAssembler):
    __doc__ = cpp.SystemAssembler.__doc__
    def __init__(self, A_form, b_form,
                 bcs=None,
                 A_coefficients=None,
                 b_coefficients=None,
                 A_function_spaces=None,
                 b_function_spaces=None,
                 cell_domains=None,
                 exterior_facet_domains=None,
                 interior_facet_domains=None,
                 form_compiler_parameters=None):
        """
        Create a SystemAssembler

        * Arguments *
           a (ufl.Form, _Form_)
              Bilinear form
           L (ufl.Form, _Form_)
              Linear form
           bcs (_DirichletBC_)
              A list or a single DirichletBC (optional)
        """
        # Extract subdomains
        subdomains = { "cell": cell_domains,
                       "exterior_facet": exterior_facet_domains,
                       "interior_facet": interior_facet_domains}

        # Wrap forms
        # First check if we got a cpp.Form which originates from cpp layer
        if isinstance(A_form, cpp.Form) and not hasattr(A_form, "_compiled_form"):

            # Then we just try to use that one
            A_dolfin_form = A_form
        else:
            A_dolfin_form = Form(A_form, A_function_spaces, A_coefficients,
                                 subdomains, form_compiler_parameters)
            
        if isinstance(b_form, cpp.Form) and not hasattr(b_form, "_compiled_form"):

            # Then we just try to use that one
            b_dolfin_form = b_form
        else:
            b_dolfin_form = Form(b_form, b_function_spaces, b_coefficients,
                                 subdomains, form_compiler_parameters)
        
        # Extract domains
        cell_domains, exterior_facet_domains, interior_facet_domains = \
                      _extract_domains(A_dolfin_form.mesh(),
                                       cell_domains,
                                       exterior_facet_domains,
                                       interior_facet_domains)

        # Attach domains to form (from Mesh, in other cases form already
        # has domain data attached)
        if cell_domains is not None:
            A_dolfin_form.set_cell_domains(cell_domains);
            b_dolfin_form.set_cell_domains(cell_domains);
        if interior_facet_domains is not None:
            A_dolfin_form.set_interior_facet_domains(interior_facet_domains);
            b_dolfin_form.set_interior_facet_domains(interior_facet_domains);
        if exterior_facet_domains is not None:
            A_dolfin_form.set_exterior_facet_domains(exterior_facet_domains);
            b_dolfin_form.set_exterior_facet_domains(exterior_facet_domains);

        # Check bcs
        bcs = _wrap_in_list(bcs, 'bcs', cpp.DirichletBC)

        # Call C++ assemble function
        cpp.SystemAssembler.__init__(self, A_dolfin_form, b_dolfin_form, bcs)
python-dolfin 1.3.0+dfsg-2 / usr / lib / python2.7 / dist-packages / dolfin / fem / assembling.py
