This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Small question on defining constant tensors

+2 votes

I currently define a tensor as follows:

m = dolfin.Constant(numpy.array([[0.98, 0, 0], [0, 0.19, 0], [0, 0, 0.19]]))

Upon inspection, this object is described as

Coefficient(TensorElement('Real', None, 0, shape=(3, 3), symmetry=None, quad_scheme=None), 54)

How do I make it clear to FEniCS that this is a diagonal (and hence symmetric) tensor?

asked Jun 9, 2016 by maartent FEniCS User (3,910 points)

1 Answer

+2 votes
 
Best answer

Hi, consider

import ufl
import dolfin.cpp as cpp
import numpy

# Following https://github.com/FEniCS/dolfin/blob/master/site-packages/dolfin/functions/constant.py

class MyTensorConstant(ufl.Coefficient, cpp.Constant):

    def __init__(self, value, cell, name=None, symmetry=None):

        cell = ufl.as_cell(cell)
        ufl_domain = None

        array = numpy.array(value)
        rank = len(array.shape)
        floats = list(map(float, array.flat))

        assert rank == 2

        ufl_element = ufl.TensorElement("Real", cell, 0, shape=array.shape,
                                        symmetry=symmetry)
        cpp.Constant.__init__(self, list(array.shape), floats)

        ufl_function_space = ufl.FunctionSpace(ufl_domain, ufl_element)
        ufl.Coefficient.__init__(self, ufl_function_space, count=self.id())

        name = name or "f_%d" % self.count()
        self.rename(name, "a Constant")

    def cell(self):
        return self.ufl_element().cell()

    def __float__(self):
        # Overriding UFL operator in this particular case.
        if self.ufl_shape:
            raise TypeError("Cannot convert nonscalar constant to float.")
        return cpp.Constant.__float__(self)

    def __str__(self):
        "x.__str__() <==> print(x)"
        return self.name()

# ----------------------------------------------------------------------------

if __name__ == '__main__':
    from dolfin import *

    mesh = UnitSquareMesh(10, 10)
    V = FunctionSpace(mesh, 'CG', 1)
    u = TrialFunction(V)
    v = TestFunction(V)

    c_values = ((1, 0), (0, 2))
    c1 = MyTensorConstant(c_values, cell=mesh.ufl_cell())
    c2 = MyTensorConstant(c_values, cell=mesh.ufl_cell(), symmetry={(0, 1): (1, 0)})

    a1 = inner(dot(c1, grad(u)), grad(v))*dx
    a2 = inner(dot(c2, grad(u)), grad(v))*dx

    print a1.signature() == a2.signature()   # Show that these are different forms

Here, a new class is introduced, following the pyDolfin's Constant class, which allows for symmetry argument to be passed to UFL's TensorConstant.

answered Jun 11, 2016 by MiroK FEniCS Expert (80,920 points)
selected Jun 26, 2016 by maartent

Thanks! Do you know by any chance how the symmetry of a UFL form is handled? More specifically, will it automatically recognize it is symmetric, or only check this at matrix level, or only invoke a solver for symmetric matrices (I use a direct solver) when a user explicitly asks for it?

Hi, with the symmetric tensor a different code is generated to assemble the matrix. The symmetry information is not used to select a direct solver - you have to do this yourself, e.g. here.

...