python-dolfin 1.3.0+dfsg-2 / usr / lib / python2.7 / dist-packages / dolfin / fem / norms.py

"""This module provides a simple way to compute various norms of
:py:class:`Vectors <dolfin.cpp.Vector>` and :py:class:`Functions
<dolfin.functions.function.Function>`, including the standard
:math:`L^2`-norm and other norms."""

# Copyright (C) 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/>.
#
# First added:  2008-03-19
# Last changed: 2009-12-11

__all__ = ["norm", "errornorm"]

from ufl import inner, grad, div, curl, dx, FiniteElement, VectorElement, Coefficient
from math import sqrt

import dolfin.cpp as cpp

from dolfin.cpp import GenericVector, GenericFunction, Function, Mesh, error, Vector
from dolfin.fem.assembling import assemble
from dolfin.fem.interpolation import interpolate
from dolfin.functions.function import Function
from dolfin.functions.functionspace import FunctionSpace, VectorFunctionSpace

dim_to_shape = {1: "interval", 2: "triangle", 3: "tetrahedron"}

def norm(v, norm_type="L2", mesh=None):
    """
    Return the norm of a given vector or function.

    *Arguments*
        v
            a :py:class:`Vector <dolfin.cpp.Vector>` or
            a :py:class:`Function <dolfin.functions.function.Function>`.
        norm_type
            see below for alternatives.
        mesh
            optional :py:class:`Mesh <dolfin.cpp.Mesh>` on
            which to compute the norm.

    If the norm type is not specified, the standard :math:`L^2` -norm
    is computed. Possible norm types include:

    *Vectors*

    ================   =================  ================
    Norm               Usage
    ================   =================  ================
    :math:`l^2`        norm(x, 'l2')      Default
    :math:`l^1`        norm(x, 'l1')
    :math:`l^\infty`   norm(x, 'linf')
    ================   =================  ================

    *Functions*

    ================  =================  =================================
    Norm              Usage              Includes the :math:`L^2` -term
    ================  =================  =================================
    :math:`L^2`       norm(v, 'L2')      Yes
    :math:`H^1`       norm(v, 'H1')      Yes
    :math:`H^1_0`     norm(v, 'H10')     No
    :math:`H` (div)   norm(v, 'Hdiv')    Yes
    :math:`H` (div)   norm(v, 'Hdiv0')   No
    :math:`H` (curl)  norm(v, 'Hcurl')   Yes
    :math:`H` (curl)  norm(v, 'Hcurl0')  No
    ================  =================  =================================

    *Examples of usage*

    .. code-block:: python

        v = Function(V)
        x = v.vector()

        print norm(x, 'linf')   # print the infinity norm of vector x

        n = norm(v)             # compute L^2 norm of v
        print norm(v, 'Hdiv')   # print H(div) norm of v
        n = norm(v, 'H1', mesh) # compute H^1 norm of v on given mesh

    """

    if not isinstance(v, (GenericVector, GenericFunction)):
        raise TypeError, "expected a GenericVector or GenericFunction"

    # Check arguments
    if not isinstance(norm_type, str):
        cpp.dolfin_error("norms.py",
                         "compute norm",
                         "Norm type must be a string, not " + str(type(norm_type)))
    if mesh is not None and not isinstance(mesh, Mesh):
        cpp.dolfin_error("norms.py",
                         "compute norm",
                         "Expecting a Mesh, not " + str(type(mesh)))

    # Select norm type
    if isinstance(v, GenericVector):
        return v.norm(norm_type.lower())

    elif (isinstance(v, Coefficient) and isinstance(v, GenericFunction)):
        if norm_type.lower() == "l2":
            M = inner(v, v)*dx()
        elif norm_type.lower() == "h1":
            M = inner(v, v)*dx() + inner(grad(v), grad(v))*dx()
        elif norm_type.lower() == "h10":
            M = inner(grad(v), grad(v))*dx()
        elif norm_type.lower() == "hdiv":
            M = inner(v, v)*dx() + div(v)*div(v)*dx()
        elif norm_type.lower() == "hdiv0":
            M = div(v)*div(v)*dx()
        elif norm_type.lower() == "hcurl":
            M = inner(v, v)*dx() + inner(curl(v), curl(v))*dx()
        elif norm_type.lower() == "hcurl0":
            M = inner(curl(v), curl(v))*dx()
        else:
            cpp.dolfin_error("norms.py",
                             "compute norm",
                             "Unknown norm type (\"%s\") for functions" % str(norm_type))
    else:
        cpp.dolfin_error("norms.py",
                         "compute norm",
                         "Unknown object type. Must be a vector or a function")

    # Get mesh
    if isinstance(v, Function) and mesh is None:
        mesh = v.function_space().mesh()

    # Assemble value
    r = assemble(M, mesh=mesh, form_compiler_parameters={"representation": "quadrature"})

    # Check value
    if r < 0.0:
        cpp.dolfin_error("norms.py",
                         "compute norm",
                         "Square of norm is negative, might be a round-off error")
    elif r == 0.0:
        return 0.0
    else:
        return sqrt(r)

def errornorm(u, uh, norm_type="l2", degree_rise=3, mesh=None):
    """
    Compute and return the error :math:`e = u - u_h` in the given norm.

    *Arguments*
        u, uh
            :py:class:`Functions <dolfin.functions.function.Function>`
        norm_type
            Type of norm. The :math:`L^2` -norm is default.
            For other norms, see :py:func:`norm <dolfin.fem.norms.norm>`.
        degree_rise
            The number of degrees above that of u_h used in the
            interpolation; i.e. the degree of piecewise polynomials used
            to approximate :math:`u` and :math:`u_h` will be the degree
            of :math:`u_h` + degree_raise.
        mesh
            Optional :py:class:`Mesh <dolfin.cpp.Mesh>` on
            which to compute the error norm.

    In simple cases, one may just define

    .. code-block:: python

        e = u - uh

    and evalute for example the square of the error in the :math:`L^2` -norm by

    .. code-block:: python

        assemble(e*e*dx(), mesh)

    However, this is not stable w.r.t. round-off errors considering that
    the form compiler may expand the expression above to::

        e*e*dx() = u*u*dx() - 2*u*uh*dx() + uh*uh*dx()

    and this might get further expanded into thousands of terms for
    higher order elements. Thus, the error will be evaluated by adding
    a large number of terms which should sum up to something close to
    zero (if the error is small).

    This module computes the error by first interpolating both
    :math:`u` and :math:`u_h` to a common space (of high accuracy),
    then subtracting the two fields (which is easy since they are
    expressed in the same basis) and then evaluating the integral.

    """

    # Check argument
    if not isinstance(u, cpp.GenericFunction):
        cpp.dolfin_error("norms.py",
                         "compute error norm",
                         "Expecting a Function or Expression for u")
    if not isinstance(uh, cpp.Function):
        cpp.dolfin_error("norms.py",
                         "compute error norm",
                         "Expecting a Function for uh")

    # Get mesh
    if isinstance(u, Function) and mesh is None:
        mesh = u.function_space().mesh()
    if isinstance(uh, Function) and mesh is None:
        mesh = uh.function_space().mesh()
    if mesh is None:
        cpp.dolfin_error("norms.py",
                         "compute error norm",
                         "Missing mesh")

    # Get rank
    if not u.value_rank() == uh.value_rank():
        cpp.dolfin_error("norms.py",
                         "compute error norm",
                         "Value ranks don't match")
    rank = u.value_rank()

    # Check that uh is associated with a finite element
    if uh.ufl_element().degree() is None:
        cpp.dolfin_error("norms.py",
                         "compute error norm",
                         "Function uh must have a finite element")

    # Degree for interpolation space. Raise degree with respect to uh.
    degree = uh.ufl_element().degree() + degree_rise

    # Check degree of 'exact' solution u
    degree_u = u.ufl_element().degree()
    if degree_u is not None and degree_u < degree:
        cpp.warning("Degree of exact solution may be inadequate for accurate result in errornorm.")

    # Create function space
    if rank == 0:
        V = FunctionSpace(mesh, "Discontinuous Lagrange", degree)
    elif rank == 1:
        V = VectorFunctionSpace(mesh, "Discontinuous Lagrange", degree)
    else:
        cpp.dolfin_error("norms.py",
                         "compute error norm",
                         "Can't handle elements of rank %d" % rank)

    # Interpolate functions into finite element space
    pi_u = interpolate(u, V)
    pi_uh = interpolate(uh, V)

    # Compute the difference
    e = Function(V)
    e.assign(pi_u)
    e.vector().axpy(-1.0, pi_uh.vector())

    # Compute norm
    return norm(e, norm_type=norm_type, mesh=mesh)