Complex Expressions with sympy operations

dolfin-adjoint
1

What is the best way to translate the process below to play nice with dolfin_adjoint?
In short, I have a sympy nsolve, and some symbolic computations.

My class ParametricGeometry generates three Expression objects that together describe a parametric surface (an extruded circular arc) and its tangent basis. This surface is parametrized by the cross sectional length l and a scalar c (relating the cross sectional width to the length). I would like to use one or both of these parameters as controls in an optimization problem.

The final Expression objects receive a bunch of arbitrarily complex formatted cpp code snippets. Ideally, I would like to rely on sympy’s symbolic differentiation as much as possible because I expect to work with more complex geometries. I.e, I know I can manually type out the final C++ code for the surface and derivatives in terms of several variables, but would like to avoid this. And this doesn’t get around the sp.nsolve.

import sympy as sp
from dolfin import Expression, pi, Constant
from dolfin_adjoint import Expression, Constant

DEGREE = 3

def ccode(z):
    return sp.printing.ccode(z)

class ParametricGeometry(object):
    def __init__(self, c=None, length=1):
        eta, L = sp.symbols('eta, length')
        eta = sp.nsolve(c**2 * eta**2 - 2*(1-sp.cos(eta)), eta, -5)

        W = c*L

        A = -0.5*W
        B = 0.5*W * (sp.sin(eta)/(1-sp.cos(eta)))
        C = 0.5*W
        D = 0.5*W * (sp.sin(eta)/(1-sp.cos(eta)))

        # center at -cL/2,
        R = L/eta
        cx = A
        cy = W*sp.sin(eta)/(2*(1-sp.cos(eta)))
        self.eta = eta
        self.radius = R.evalf(subs={L: 1})
        self.center = [cx.evalf(subs={L: length}), cy.evalf(subs={L: length})]

        x1, x2 = sp.symbols('x[0], x[1]')
        gamma_sp = [None, None, None]
        gamma_sp[0] = A*sp.cos(x1*L/R) + B*sp.sin(x1*L/R) + C
        gamma_sp[1] = x2
        gamma_sp[2] = A*sp.sin(x1*L/R) - B*sp.cos(x1*L/R) + D
        self._build_geo(gamma_sp, length, x1, x2)

    def _build_geo(self, gamma_sp, length, x1, x2):

        # Construct the actual fenics map X = gamma(xi)
        self.gamma = Expression((ccode(gamma_sp[0]).replace('M_PI', 'pi'),
                                 ccode(gamma_sp[1]).replace('M_PI', 'pi'),
                                 ccode(gamma_sp[2]).replace('M_PI', 'pi')),
                                pi=pi,
                                length=length,
                                degree=DEGREE,
                                name="gamma")

        # Get the induced base from the mapping to the referential geometry
        # G_i = ∂X/xi_i^2
        self.Gsub1 = Expression((ccode(gamma_sp[0].diff(x1)).replace('M_PI', 'pi'),
                                 ccode(gamma_sp[1].diff(x1)).replace('M_PI', 'pi'),
                                 ccode(gamma_sp[2].diff(x1)).replace('M_PI', 'pi')),
                                pi=pi,
                                length=length,
                                degree=DEGREE,
                                name="Gsub1")

        self.Gsub2 = Expression((ccode(gamma_sp[0].diff(x2)).replace('M_PI', 'pi'),
                                 ccode(gamma_sp[1].diff(x2)).replace('M_PI', 'pi'),
                                 ccode(gamma_sp[2].diff(x2)).replace('M_PI', 'pi')),
                                pi=pi,
                                length=length,
                                degree=DEGREE,
                                name="Gsub2")