Source code for simudo.fem.ufl_extra

from __future__ import absolute_import, division, print_function

import math
import unittest
from builtins import bytes, dict, int, range, str, super

import dolfin
import ufl
from ufl.constantvalue import (
    FloatValue, IntValue, ScalarValue, Zero, as_ufl, is_true_ufl_scalar)
from ufl.core.ufl_type import ufl_type
from ufl.mathfunctions import MathFunction

__all__ = ['ExpM1', 'expm1']

'''
note: MathFunction.derivative() is an undocumented attribute which can
be used to specify a function's derivative; see
`ufl.algorithms.apply_derivatives.GenericDerivativeRuleset.math_function`
for usage.
'''

def ufl_mathfunction(f, cls):
    # taken from ufl.operators
    f = as_ufl(f)
    r = cls(f)
    if isinstance(r, (ScalarValue, Zero, int, float)):
        return float(r)
    return r

[docs]@ufl_type() class ExpM1(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (ScalarValue, Zero)): return FloatValue(math.expm1(float(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "expm1", argument)
[docs] def derivative(self): f, = self.ufl_operands return ufl.exp(f)
[docs]def expm1(x): return ufl_mathfunction(x, ExpM1)
class Test(unittest.TestCase): def setUp(self): self.mesh = mesh = dolfin.UnitIntervalMesh(30) self.element = element = dolfin.FiniteElement("CG", mesh.ufl_cell(), 3) self.W = W = dolfin.FunctionSpace(mesh, element) self.u = dolfin.Function(W, name='u') self.v = dolfin.TestFunction(W) self.x = dolfin.SpatialCoordinate(mesh)[0] def test_expm1(self): x = self.x u1 = dolfin.project(expm1(x), self.W) u2 = dolfin.project(dolfin.exp(x)-1.0) err = dolfin.assemble((u1-u2)**2*dolfin.dx) self.assertLessEqual(err, 1e-8) def test_expm1_derivative(self): ''' check that we can take functional derivatives ''' x = self.x u = self.u v = self.v dolfin.project(x, self.W, function=u) F = expm1(u*u)*v*dolfin.dx - (dolfin.exp(x) - 1.0)*v*dolfin.dx from ..newtonsolver import NewtonSolver solver = NewtonSolver(F, u, []) solver.parameters['relative_tolerance'] = 1e-14 solver.solve() u_e = dolfin.sqrt(x) err = dolfin.assemble((u - u_e)**2*dolfin.dx) self.assertLessEqual(err, 1e-8)