Source code for sympy.galgebra.manifold

# sympy/galgebra/manifold.py

"""
manifold.py defines the Manifold class which allows one to create a
vector manifold (manifold defined by vector field of coordinates in
embedding vector space) calculate the tangent vectors and derivatives
of tangent vectors.

Once manifold is created multivector fields can be constructed in the
tangent space and all the geometric algebra products and derivatives
of the multivector fields calculated.

Note that all calculations are done in the embedding space.  Future
versions of the code will allow manifolds defined purely in terms of
a metric.
"""

from __future__ import print_function

from itertools import combinations
from os import system
import copy

from sympy import trigsimp, simplify

from sympy.galgebra.ga import MV
from sympy.galgebra.debug import oprint
from sympy.galgebra.ncutil import linear_expand
from sympy.galgebra.printing import find_executable


def fct_to_str(fct_names):
    import sys
    current_file = open(sys.argv[0], 'r')
    file_str = current_file.read()
    current_file.close()

    if isinstance(fct_names, str):
        return fct_names

    fcts_str = ''

    for fct_name in fct_names:
        start_def = file_str.find('\ndef ' + fct_name)
        end_def = file_str.find('\ndef ', start_def + 5)
        start_class = file_str.find('\nclass ', start_def + 5)
        end_def = min(end_def, start_class)
        fcts_str += file_str[start_def:end_def]
    return fcts_str


def VectorComponents(X, basis):
    (coefs, bases) = linear_expand(X.obj)
    cdict = {}
    for (coef, base) in zip(coefs, bases):
        cdict[str(base)] = coef
    comp = []
    for base in basis:
        if base in cdict:
            comp.append(cdict[base])
        else:
            comp.append(0)
    return comp


def FillTemplate(self, template):
    Nd = 0
    var = []
    id_old = 0
    while True:
        id_new = template.find('$', id_old + 1)
        if id_new == -1:
            break
        Nd += 1
        if Nd % 2 == 0:
            var.append(template[id_old + 1:id_new])
        id_old = id_new

    var.sort(reverse=True)

    for v in var:
        template = template.replace('$' + v + '$', str(eval('self.' + v)))

    return template


[docs]class Manifold: def __init__(self, x, coords, debug=False, I=None): """ coords: list of coordinate variables x: vector fuction of coordinate variables (parametric surface) """ self.I = I self.x = x self.coords = coords self.basis = [] self.basis_str = [] self.embedded_basis = [] for u in coords: tv = x.diff(u) self.basis.append(tv) (coefs, bases) = linear_expand(tv.obj) tc = {} for (coef, base) in zip(coefs, bases): str_base = str(base) tc[str_base] = coef if str_base not in self.embedded_basis: self.embedded_basis.append(str_base) self.basis_str.append(tc) self.gij = [] for base1 in self.basis: tmp = [] for base2 in self.basis: tmp.append(simplify(trigsimp((base1 | base2).scalar()))) self.gij.append(tmp) for tv in self.basis_str: for base in self.embedded_basis: if base not in tv: tv[base] = 0 self.dim = len(self.basis) indexes = tuple(range(self.dim)) self.index = [()] for i in indexes: self.index.append(tuple(combinations(indexes, i + 1))) self.index = tuple(self.index) self.MFbasis = [[MV.ONE], self.basis] for igrade in self.index[2:]: grade = [] for iblade in igrade: blade = MV(1, 'scalar') for ibasis in iblade: blade ^= self.basis[ibasis] blade = blade.trigsimp(deep=True, recursive=True) grade.append(blade) self.MFbasis.append(grade) self.E = self.MFbasis[-1][0] self.E_sq = trigsimp((self.E * self.E).scalar(), deep=True, recursive=True) duals = copy.copy(self.MFbasis[-2]) duals.reverse() sgn = 1 self.rbasis = [] for dual in duals: recpv = (sgn * dual * self.E).trigsimp(deep=True, recursive=True) self.rbasis.append(recpv) sgn = -sgn self.dbasis = [] for base in self.basis: dbase = [] for coord in self.coords: d = base.diff(coord).trigsimp(deep=True, recursive=True) dbase.append(d) self.dbasis.append(dbase) self.surface = {} (coefs, bases) = linear_expand(self.x.obj) for (coef, base) in zip(coefs, bases): self.surface[str(base)] = coef self.grad = MV() self.grad.is_grad = True self.grad.blade_rep = True self.grad.igrade = 1 self.grad.rcpr_bases_MV = [] for rbase in self.rbasis: self.grad.rcpr_bases_MV.append(rbase / self.E_sq) self.grad.rcpr_bases_MV = tuple(self.grad.rcpr_bases_MV) self.grad.coords = self.coords self.grad.norm = self.E_sq self.grad.connection = {} if debug: oprint('x', self.x, 'coords', self.coords, 'basis vectors', self.basis, 'index', self.index, 'basis blades', self.MFbasis, 'E', self.E, 'E**2', self.E_sq, '*basis', duals, 'rbasis', self.rbasis, 'basis derivatives', self.dbasis, 'surface', self.surface, 'basis strings', self.basis_str, 'embedding basis', self.embedded_basis, 'metric tensor', self.gij) def Basis(self): return tuple(self.basis) def Grad(self, F): # Intrisic Derivative dF = 0 for (rbase, coord) in zip(self.rbasis, self.coords): dF += rbase * F.diff(coord) dF = dF.simplify() dF = dF / self.E_sq return dF def D(self, F): # Covariant Derivative dF = self.Grad(F) return self.Proj(dF) def S(self, a): # Shape Tensor return def Proj(self, F): PF = (F < self.E) * self.E PF = PF.simplify() PF = PF.trigsimp(deep=True, recursive=True) return (PF / self.E_sq).simplify() def Reject(self, F): return (F - self.Proj(F)).simplify() def DD(self, v, f, opstr=False): mf_comp = [] for e in self.rbasis: mf_comp.append((v | e).scalar() / self.E_sq) result = MV() op = '' for (coord, comp) in zip(self.coords, mf_comp): result += comp * (f.diff(coord)) if opstr: op += '(' + str(comp) + ')D{' + str(coord) + '}+' if opstr: return str(result), op[:-1] return result def Plot2DSurface(self, u_range, v_range, surf=True, grid=True, tan=1.0, scalar_field=None, skip=[1, 1], fct_def=None): plot_template = \ """ from numpy import mgrid,shape,swapaxes,zeros,log,exp,sin,cos,tan $fct_def$ eps = 1.0e-6 u_r = $u_range$ v_r = $v_range$ $coords$ = mgrid[u_r[0]:u_r[1]+eps:(u_r[1]-u_r[0])/float(u_r[2]-1),\\ v_r[0]:v_r[1]+eps:(v_r[1]-v_r[0])/float(v_r[2]-1)] X = $surface$ scal_tan = $tan$ x = X['ex'] y = X['ey'] z = X['ez'] du = $basis_str[0]$ dv = $basis_str[1]$ Zero = zeros(shape(x)) if scal_tan > 0.0: du_x = Zero+du['ex'] du_y = Zero+du['ey'] du_z = Zero+du['ez'] dv_x = Zero+dv['ex'] dv_y = Zero+dv['ey'] dv_z = Zero+dv['ez'] f = $scalar_field$ n = $n$ skip = $skip$ su = skip[0] sv = skip[1] if f[0] != None: dn_x = f[0]*n[0] dn_y = f[0]*n[1] dn_z = f[0]*n[2] from mayavi.mlab import plot3d,quiver3d,mesh,figure figure(bgcolor=(1.0,1.0,1.0)) if $surf$: mesh(x,y,z,colormap="gist_earth") if $grid$: for i in range(shape(u)[0]): plot3d(x[i,],y[i,],z[i,],line_width=1.0,color=(0.0,0.0,0.0),tube_radius=None) xr = swapaxes(x,0,1) yr = swapaxes(y,0,1) zr = swapaxes(z,0,1) for i in range(shape(u)[1]): plot3d(xr[i,],yr[i,],zr[i,],line_width=1.0,color=(0.0,0.0,0.0),tube_radius=None) if scal_tan > 0.0: quiver3d(x[::su,::sv],y[::su,::sv],z[::su,::sv],\\ du_x[::su,::sv],du_y[::su,::sv],du_z[::su,::sv],scale_factor=scal_tan,\\ line_width=1.0,color=(0.0,0.0,0.0),scale_mode='vector',mode='arrow',resolution=16) quiver3d(x[::su,::sv],y[::su,::sv],z[::su,::sv],\\ dv_x[::su,::sv],dv_y[::su,::sv],dv_z[::su,::sv],scale_factor=scal_tan,\\ line_width=1.0,color=(0.0,0.0,0.0),scale_mode='vector',mode='arrow',resolution=16) if f[0] != None: quiver3d(x[::su,::sv],y[::su,::sv],z[::su,::sv],\\ dn_x[::su,::sv],dn_y[::su,::sv],dn_z[::su,::sv],\\ line_width=1.0,color=(0.0,0.0,0.0),scale_mode='none',mode='cone',\\ resolution=16,opacity=0.5) """ if len(self.coords) != 2: return self.skip = skip self.surf = surf self.grid = grid self.tan = tan if fct_def is None: self.fct_def = ' ' else: self.fct_def = fct_to_str(fct_def) self.u_range = u_range self.v_range = v_range self.scalar_field = [scalar_field] print(self.I, '\n', self.basis[0], '\n', self.basis[1]) self.normal = -self.I * (self.basis[0] ^ self.basis[1]) self.n = VectorComponents(self.normal, ['ex', 'ey', 'ez']) msurf = open('manifold_surf.py', 'w') plot_template = FillTemplate(self, plot_template) msurf.write(plot_template) msurf.close() mayavi2 = find_executable('mayavi2') if mayavi2 is None: return system(mayavi2 + ' manifold_surf.py &') return