gradient_flux.py 4.95 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
1 2 3 4 5 6 7 8 9 10 11 12 13
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Feb 24 14:36:40 2020

@author: bleyerj
"""
from dolfin import *
import numpy as np
import ufl
from mfront_wrapper.utils import *

class Gradient:
14
    def __init__(self, variable, expression, name=None, symmetric=None):
Jeremy BLEYER's avatar
Jeremy BLEYER committed
15
        self.variable = variable
16 17 18
        if symmetric is None:
            self.expression = expression
        elif symmetric:
19 20 21 22
            if self.variable.geometric_dimension()==2:
                self.expression = as_vector([symmetric_tensor_to_vector(expression)[i] for i in range(4)])
            else:
                self.expression = symmetric_tensor_to_vector(expression)
23
        else:
24 25 26 27
            if self.variable.geometric_dimension()==2:
                self.expression = as_vector([nonsymmetric_tensor_to_vector(expression)[i] for i in range(5)])
            else:
                self.expression = nonsymmetric_tensor_to_vector(expression)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
28 29 30 31
        shape = ufl.shape(self.expression)
        if len(shape)==1:
            self.shape = shape[0]
        elif shape==():
32
            self.shape = 0
33
            #self.expression = as_vector([self.expression])
Jeremy BLEYER's avatar
Jeremy BLEYER committed
34 35 36 37 38 39 40 41 42 43 44 45
        else:
            self.shape = shape
        if name is None:
            self.name = "Gradient-{}".format(np.random.randint(1e6))
        else:
            self.name = name
    def __call__(self, v):
        return ufl.replace(self.expression, {self.variable:v})
    def variation(self, v):
        """ Directional derivative in direction v """
        return ufl.derivative(self.expression, self.variable, v)

46 47 48 49 50 51 52 53
    def initialize_function(self, mesh, quadrature_degree):
        We = get_quadrature_element(mesh.ufl_cell(), quadrature_degree, self.shape)
        self.function_space = FunctionSpace(mesh, We)
        self.function = Function(self.function_space, name=self.name)

    def update(self, x):
        self.function.vector().set_local(x)

Jeremy BLEYER's avatar
Jeremy BLEYER committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
class Grad(Gradient):
    """ Gradient operator : grad(u) """
    def __init__(self, variable, name=None):
        Gradient.__init__(self, variable,
                          nonsymmetric_tensor_to_vector(grad(variable)), name)

class GradT(Gradient):
    """ Transposed gradient operator : grad(u).T """
    def __init__(self, variable, name=None):
        Gradient.__init__(self, variable,
                          nonsymmetric_tensor_to_vector(grad(variable).T), name)

class SymGrad(Gradient):
    """ Symmetric gradient operator : sym(grad(u)) """
    def __init__(self, variable, name=None):
        Gradient.__init__(self, variable,
                          symmetric_tensor_to_vector(sym(grad(variable))), name)

class Grad_Id(Gradient):
    """ Transformation gradient operator : grad(u) + Id(d) """
    def __init__(self, variable, name=None):
        Gradient.__init__(self, variable,
                          symmetric_tensor_to_vector(grad(variable) +
                          Identity(variable.geometric_dimension()), T22=1), name)
class Var(Gradient):
    """ The variable itself : u """
    def __init__(self, variable, name=None):
        Gradient.__init__(self, variable, variable, name)


class Flux:
85 86 87 88
    def __init__(self, name, shape, gradients=[]):
        self.shape = shape
        self.name = name
        self.gradients = gradients
Jeremy BLEYER's avatar
Jeremy BLEYER committed
89 90 91 92 93

    def initialize_functions(self, mesh, quadrature_degree):
        We = get_quadrature_element(mesh.ufl_cell(), quadrature_degree, self.shape)
        W = FunctionSpace(mesh, We)
        self.function = Function(W, name=self.name)
94
        values = [Function(FunctionSpace(mesh,
Jeremy BLEYER's avatar
Jeremy BLEYER committed
95
                          get_quadrature_element(mesh.ufl_cell(),
96 97
                          quadrature_degree, dim=(self.shape, g.shape))),
                           name="d{}_d{}".format(self.name, g.name))
Jeremy BLEYER's avatar
Jeremy BLEYER committed
98
                          for g in self.gradients]
99 100
        keys = [g.name for g in self.gradients]
        self.tangent_blocks = dict(zip(keys, values))
101 102 103 104

    def update(self, x):
        self.function.vector().set_local(x)

105 106 107 108 109 110
    def previous(self):
        try:
            self.previous = Function(self.function.function_space(), name=self.name+"_previous")
            self.previous.assign(self.function)
            return self.previous
        except:
111
            raise ValueError("Function must be initialized first.")
Jeremy BLEYER's avatar
Jeremy BLEYER committed
112

113 114 115
# class Stress(Flux):
#     def __init__(self, variable, name="sigma"):
#         Flux.__init__(self, SymGrad(variable, name="eps"), name=name)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
116 117


118 119 120 121
# class ThermalFlux(Flux):
#     def __init__(self, variable, name="j"):
#         Flux.__init__(self, [Grad(variable, name="grad(T)"), Var(variable, name="T")],
#                              variable.geometric_dimension(), name=name)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
122

123 124 125 126 127 128 129 130 131 132 133 134 135
# mesh = UnitSquareMesh(10, 10)
# V = VectorFunctionSpace(mesh, "CG", 1)
# VT = FunctionSpace(mesh, "CG", 1)

# u = Function(V, name="u")
# T = Function(VT, name="Temperature")
# v = Function(V, name="v")
# g = Grad_Id(u, name="F")
# sig = Stress(u)
# j = ThermalFlux(T)
# j.initialize_functions(mesh, 2)
# print([ufl.shape(tang) for tang in j.tangents])
# print([tang.name() for tang in j.tangents])