mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-09 01:35:52 +08:00
Merge branch 'develop' of https://github.com/simpeg/simpeg into cylClean
This commit is contained in:
rowanc1
@@ -239,7 +239,7 @@ class BaseMesh(object):
|
||||
:return: projected face vector
|
||||
"""
|
||||
assert type(fV) == np.ndarray, 'fV must be an ndarray'
|
||||
assert len(fV.shape) == 2 and fV.shape[0] == np.sum(self.nF) and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
|
||||
assert len(fV.shape) == 2 and fV.shape[0] == self.nF and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
|
||||
return np.sum(fV*self.normals, 1)
|
||||
|
||||
def projectEdgeVector(self, eV):
|
||||
@@ -251,7 +251,7 @@ class BaseMesh(object):
|
||||
:return: projected edge vector
|
||||
"""
|
||||
assert type(eV) == np.ndarray, 'eV must be an ndarray'
|
||||
assert len(eV.shape) == 2 and eV.shape[0] == np.sum(self.nE) and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
|
||||
assert len(eV.shape) == 2 and eV.shape[0] == self.nE and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
|
||||
return np.sum(eV*self.tangents, 1)
|
||||
|
||||
|
||||
|
||||
@@ -142,7 +142,14 @@ class InnerProducts(object):
|
||||
Note that this is completed for each cell in the mesh at the same time.
|
||||
|
||||
"""
|
||||
if M.dim == 2:
|
||||
if M.dim == 1:
|
||||
v = np.sqrt(0.5*M.vol)
|
||||
V1 = sdiag(v) # We will multiply on each side to keep symmetry
|
||||
|
||||
Px = _getFacePx(M)
|
||||
P000 = V1*Px('fXm')
|
||||
P100 = V1*Px('fXp')
|
||||
elif M.dim == 2:
|
||||
# Square root of cell volume multiplied by 1/4
|
||||
v = np.sqrt(0.25*M.vol)
|
||||
V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry
|
||||
@@ -168,9 +175,13 @@ class InnerProducts(object):
|
||||
P111 = V3*Pxxx('fXp', 'fYp', 'fZp')
|
||||
|
||||
Mu = _makeTensor(M, mu)
|
||||
A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110
|
||||
P = [P000, P100, P010, P110]
|
||||
if M.dim == 3:
|
||||
A = P000.T*Mu*P000 + P100.T*Mu*P100
|
||||
P = [P000, P100]
|
||||
|
||||
if M.dim > 1:
|
||||
A = A + P010.T*Mu*P010 + P110.T*Mu*P110
|
||||
P += [P010, P110]
|
||||
if M.dim > 2:
|
||||
A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
|
||||
P += [P001, P101, P011, P111]
|
||||
if returnP:
|
||||
@@ -246,8 +257,10 @@ class InnerProducts(object):
|
||||
Note that this is completed for each cell in the mesh at the same time.
|
||||
|
||||
"""
|
||||
if M.dim == 1:
|
||||
raise NotImplementedError('getEdgeInnerProduct not implemented for 1D')
|
||||
# We will multiply by V on each side to keep symmetry
|
||||
if M.dim == 2:
|
||||
elif M.dim == 2:
|
||||
# Square root of cell volume multiplied by 1/4
|
||||
v = np.sqrt(0.25*M.vol)
|
||||
V = sdiag(np.r_[v, v])
|
||||
@@ -306,7 +319,13 @@ def _makeTensor(M, sigma):
|
||||
elif type(sigma) is float:
|
||||
sigma = np.ones(self.nC)*sigma
|
||||
|
||||
if M.dim == 2:
|
||||
if M.dim == 1:
|
||||
if sigma.size == M.nC: # Isotropic!
|
||||
sigma = mkvc(sigma) # ensure it is a vector.
|
||||
Sigma = sdiag(sigma)
|
||||
else:
|
||||
raise Exception('Unexpected shape of sigma')
|
||||
elif M.dim == 2:
|
||||
if sigma.size == M.nC: # Isotropic!
|
||||
sigma = mkvc(sigma) # ensure it is a vector.
|
||||
Sigma = sdiag(np.r_[sigma, sigma])
|
||||
@@ -316,6 +335,8 @@ def _makeTensor(M, sigma):
|
||||
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
|
||||
row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
|
||||
Sigma = sp.vstack((row1, row2))
|
||||
else:
|
||||
raise Exception('Unexpected shape of sigma')
|
||||
elif M.dim == 3:
|
||||
if sigma.size == M.nC: # Isotropic!
|
||||
sigma = mkvc(sigma) # ensure it is a vector.
|
||||
@@ -327,8 +348,15 @@ def _makeTensor(M, sigma):
|
||||
row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
|
||||
row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
|
||||
Sigma = sp.vstack((row1, row2, row3))
|
||||
else:
|
||||
raise Exception('Unexpected shape of sigma')
|
||||
return Sigma
|
||||
|
||||
|
||||
def _getFacePx(M):
|
||||
assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
|
||||
return _getFacePx_Rectangular(M)
|
||||
|
||||
def _getFacePxx(M):
|
||||
if M._meshType == 'TREE':
|
||||
return M._getFacePxx
|
||||
@@ -353,6 +381,23 @@ def _getEdgePxxx(M):
|
||||
|
||||
return _getEdgePxxx_Rectangular(M)
|
||||
|
||||
def _getFacePx_Rectangular(M):
|
||||
"""Returns a function for creating projection matrices
|
||||
|
||||
"""
|
||||
ii = np.int64(range(M.nCx))
|
||||
|
||||
def Px(xFace):
|
||||
"""
|
||||
xFace is 'fXp' or 'fXm'
|
||||
"""
|
||||
posFx = 0 if xFace == 'fXm' else 1
|
||||
IND = ii + posFx
|
||||
PX = sp.csr_matrix((np.ones(M.nC), (range(M.nC), IND)), shape=(M.nC, M.nF))
|
||||
return PX
|
||||
|
||||
return Px
|
||||
|
||||
def _getFacePxx_Rectangular(M):
|
||||
"""returns a function for creating projection matrices
|
||||
|
||||
@@ -410,7 +455,7 @@ def _getFacePxx_Rectangular(M):
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
|
||||
PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nF)))
|
||||
PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
|
||||
|
||||
if M._meshType == 'LOM':
|
||||
I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
|
||||
@@ -467,7 +512,7 @@ def _getFacePxxx_Rectangular(M):
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nF))).tocsr()
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
|
||||
|
||||
if M._meshType == 'LOM':
|
||||
I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
|
||||
@@ -502,7 +547,7 @@ def _getEdgePxx_Rectangular(M):
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
|
||||
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nE))).tocsr()
|
||||
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
|
||||
|
||||
if M._meshType == 'LOM':
|
||||
I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
|
||||
@@ -545,7 +590,7 @@ def _getEdgePxxx_Rectangular(M):
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nE))).tocsr()
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
|
||||
|
||||
if M._meshType == 'LOM':
|
||||
I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
|
||||
|
||||
@@ -32,6 +32,7 @@ class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
|
||||
|
||||
def __init__(self, nodes):
|
||||
assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
|
||||
assert len(nodes) > 1, "len(node) must be greater than 1"
|
||||
|
||||
for i, nodes_i in enumerate(nodes):
|
||||
assert type(nodes_i) == np.ndarray, ("nodes[%i] is not a numpy array." % i)
|
||||
|
||||
@@ -679,6 +679,8 @@ class TreeMesh(InnerProducts, BaseMesh):
|
||||
|
||||
def __init__(self, h_in, x0=None):
|
||||
assert type(h_in) is list, 'h_in must be a list'
|
||||
assert len(h_in) > 1, "len(h_in) must be greater than 1"
|
||||
|
||||
h = range(len(h_in))
|
||||
for i, h_i in enumerate(h_in):
|
||||
if type(h_i) in [int, long, float]:
|
||||
|
||||
@@ -238,6 +238,100 @@ class TensorView(object):
|
||||
if showIt: plt.show()
|
||||
return ph
|
||||
|
||||
def plotSlice(self, v, vType='CC',
|
||||
normal='Z', ind=None, grid=False, view='real',
|
||||
ax=None, clim=None, showIt=False,
|
||||
gridOpts={'color':'k'}
|
||||
):
|
||||
viewOpts = ['real','imag','abs','vec','|vec|']
|
||||
normalOpts = ['X', 'Y', 'Z']
|
||||
vTypeOpts = ['CC','F','E']
|
||||
|
||||
# Some user error checking
|
||||
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
|
||||
assert self.dim == 3, 'Must be a 3D mesh.'
|
||||
assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
|
||||
assert normal in normalOpts, "normal must be in ['%s']" % "','".join(normalOpts)
|
||||
assert type(grid) is bool, 'grid must be a boolean'
|
||||
|
||||
szSliceDim = getattr(self, 'nC'+normal.lower()) #: Size of the sliced dimension
|
||||
if ind is None: ind = int(szSliceDim/2)
|
||||
assert type(ind) in [int, long], 'ind must be an integer'
|
||||
|
||||
normalizeVector = False
|
||||
if view == '|vec|':
|
||||
view = 'vec'
|
||||
normalizeVector = True
|
||||
|
||||
if ax is None:
|
||||
fig = plt.figure(1)
|
||||
fig.clf()
|
||||
ax = plt.subplot(111)
|
||||
else:
|
||||
assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes"
|
||||
fig = ax.figure
|
||||
|
||||
# The slicing and plotting code!!
|
||||
|
||||
def getIndSlice(v):
|
||||
if normal == 'X': v = v[ind,:,:]
|
||||
elif normal == 'Y': v = v[:,ind,:]
|
||||
elif normal == 'Z': v = v[:,:,ind]
|
||||
return v
|
||||
|
||||
def doSlice(v):
|
||||
if vType == 'CC':
|
||||
return getIndSlice(self.r(v,'CC','CC','M'))
|
||||
# Now just deal with 'F' and 'E'
|
||||
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
|
||||
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
|
||||
if view == 'vec':
|
||||
v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M')
|
||||
outSlice = []
|
||||
if 'X' not in normal: outSlice.append(getIndSlice(v[0]))
|
||||
if 'Y' not in normal: outSlice.append(getIndSlice(v[1]))
|
||||
if 'Z' not in normal: outSlice.append(getIndSlice(v[2]))
|
||||
return outSlice
|
||||
else:
|
||||
return getIndSlice(self.r(v,'CC','CC','M'))
|
||||
|
||||
h2d = []
|
||||
if 'X' not in normal: h2d.append(self.hx)
|
||||
if 'Y' not in normal: h2d.append(self.hy)
|
||||
if 'Z' not in normal: h2d.append(self.hz)
|
||||
tM = self.__class__(h2d) #: Temp Mesh
|
||||
|
||||
out = ()
|
||||
if view in ['real','imag','abs']:
|
||||
v = getattr(np,view)(v) # e.g. np.real(v)
|
||||
v = doSlice(v)
|
||||
if clim is None:
|
||||
clim = [v.min(),v.max()]
|
||||
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1]),)
|
||||
elif view in ['vec']:
|
||||
U, V = doSlice(v)
|
||||
if clim is None:
|
||||
clim = [v.min(),v.max()]
|
||||
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, 0.5*(U+V).T, vmin=clim[0], vmax=clim[1]),)
|
||||
if normalizeVector:
|
||||
length = np.sqrt(U**2+V**2)
|
||||
U, V = U/length, V/length
|
||||
out += (quiver( tM.gridCC[:,0], tM.gridCC[:,1], U, V),)
|
||||
|
||||
if grid:
|
||||
xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten()
|
||||
xYGrid = np.c_[tM.vectorNy[0]*np.ones(tM.nNx),tM.vectorNy[-1]*np.ones(tM.nNx),np.nan*np.ones(tM.nNx)].flatten()
|
||||
yXGrid = np.c_[tM.vectorNx[0]*np.ones(tM.nNy),tM.vectorNx[-1]*np.ones(tM.nNy),np.nan*np.ones(tM.nNy)].flatten()
|
||||
yYGrid = np.c_[tM.vectorNy,tM.vectorNy,np.nan*np.ones(tM.nNy)].flatten()
|
||||
out += (ax.plot(np.r_[xXGrid,yXGrid],np.r_[xYGrid,yYGrid],**gridOpts)[0],)
|
||||
|
||||
ax.set_xlabel('y' if normal == 'X' else 'x')
|
||||
ax.set_ylabel('y' if normal == 'Z' else 'z')
|
||||
ax.set_title('Slice %d' % ind)
|
||||
|
||||
if showIt: plt.show()
|
||||
return out
|
||||
|
||||
def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False):
|
||||
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
|
||||
|
||||
|
||||
@@ -3,32 +3,6 @@ import unittest
|
||||
from TestUtils import OrderTest
|
||||
|
||||
|
||||
# MATLAB code:
|
||||
|
||||
# syms x y z
|
||||
|
||||
# ex = x.^2+y.*z;
|
||||
# ey = (z.^2).*x+y.*z;
|
||||
# ez = y.^2+x.*z;
|
||||
|
||||
# e = [ex;ey;ez];
|
||||
|
||||
# sigma1 = x.*y+1;
|
||||
# sigma2 = x.*z+2;
|
||||
# sigma3 = 3+z.*y;
|
||||
# sigma4 = 0.1.*x.*y.*z;
|
||||
# sigma5 = 0.2.*x.*y;
|
||||
# sigma6 = 0.1.*z;
|
||||
|
||||
# S1 = [sigma1,0,0;0,sigma1,0;0,0,sigma1];
|
||||
# S2 = [sigma1,0,0;0,sigma2,0;0,0,sigma3];
|
||||
# S3 = [sigma1,sigma4,sigma5;sigma4,sigma2,sigma6;sigma5,sigma6,sigma3];
|
||||
|
||||
# i1 = int(int(int(e.'*S1*e,x,0,1),y,0,1),z,0,1);
|
||||
# i2 = int(int(int(e.'*S2*e,x,0,1),y,0,1),z,0,1);
|
||||
# i3 = int(int(int(e.'*S3*e,x,0,1),y,0,1),z,0,1);
|
||||
|
||||
|
||||
class TestInnerProducts(OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
@@ -54,14 +28,14 @@ class TestInnerProducts(OrderTest):
|
||||
Gc = self.M.gridCC
|
||||
if self.sigmaTest == 1:
|
||||
sigma = np.c_[call(sigma1, Gc)]
|
||||
analytic = 647./360 # Found using matlab symbolic toolbox.
|
||||
analytic = 647./360 # Found using sympy.
|
||||
elif self.sigmaTest == 3:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
analytic = 37./12 # Found using matlab symbolic toolbox.
|
||||
analytic = 37./12 # Found using sympy.
|
||||
elif self.sigmaTest == 6:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc),
|
||||
call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)]
|
||||
analytic = 69881./21600 # Found using matlab symbolic toolbox.
|
||||
analytic = 69881./21600 # Found using sympy.
|
||||
|
||||
if self.location == 'edges':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
|
||||
@@ -143,13 +117,13 @@ class TestInnerProducts2D(OrderTest):
|
||||
Gc = self.M.gridCC
|
||||
if self.sigmaTest == 1:
|
||||
sigma = np.c_[call(sigma1, Gc)]
|
||||
analytic = 144877./360 # Found using matlab symbolic toolbox. z=5
|
||||
analytic = 144877./360 # Found using sympy. z=5
|
||||
elif self.sigmaTest == 2:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
|
||||
analytic = 189959./120 # Found using matlab symbolic toolbox. z=5
|
||||
analytic = 189959./120 # Found using sympy. z=5
|
||||
elif self.sigmaTest == 3:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
analytic = 781427./360 # Found using matlab symbolic toolbox. z=5
|
||||
analytic = 781427./360 # Found using sympy. z=5
|
||||
|
||||
if self.location == 'edges':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
|
||||
@@ -206,5 +180,97 @@ class TestInnerProducts2D(OrderTest):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
|
||||
class TestInnerProducts1D(OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 1
|
||||
meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
|
||||
y = 12 # Because 12 is just such a great number.
|
||||
z = 5 # Because 5 is just such a great number as well!
|
||||
|
||||
call = lambda fun, x: fun(x)
|
||||
|
||||
ex = lambda x: x**2+y*z
|
||||
|
||||
sigma1 = lambda x: x*y+1
|
||||
|
||||
Gc = self.M.gridCC
|
||||
sigma = call(sigma1, Gc)
|
||||
analytic = 128011./5 # Found using sympy. y=12, z=5
|
||||
|
||||
if self.location == 'faces':
|
||||
F = call(ex, self.M.gridFx)
|
||||
A = self.M.getFaceInnerProduct(sigma)
|
||||
numeric = F.T.dot(A.dot(F))
|
||||
|
||||
err = np.abs(numeric - analytic)
|
||||
return err
|
||||
|
||||
def test_order1_faces(self):
|
||||
self.name = "1D Face Inner Product"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.orderTest()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
if __name__ == '__main__' and False:
|
||||
import sympy
|
||||
|
||||
x,y,z = sympy.symbols(['x','y','z'])
|
||||
ex = x**2+y*z
|
||||
ey = (z**2)*x+y*z
|
||||
ez = y**2+x*z
|
||||
e = sympy.Matrix([ex,ey,ez])
|
||||
|
||||
sigma1 = x*y+1
|
||||
sigma2 = x*z+2
|
||||
sigma3 = 3+z*y
|
||||
sigma4 = 0.1*x*y*z
|
||||
sigma5 = 0.2*x*y
|
||||
sigma6 = 0.1*z
|
||||
|
||||
S1 = sympy.Matrix([[sigma1,0,0],[0,sigma1,0],[0,0,sigma1]])
|
||||
S2 = sympy.Matrix([[sigma1,0,0],[0,sigma2,0],[0,0,sigma3]])
|
||||
S3 = sympy.Matrix([[sigma1,sigma4,sigma5],[sigma4,sigma2,sigma6],[sigma5,sigma6,sigma3]])
|
||||
|
||||
print '3D'
|
||||
print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
|
||||
|
||||
z = 5
|
||||
ex = x**2+y*z
|
||||
ey = (z**2)*x+y*z
|
||||
e = sympy.Matrix([ex,ey])
|
||||
|
||||
sigma1 = x*y+1
|
||||
sigma2 = x*z+2
|
||||
sigma3 = 3+z*y
|
||||
|
||||
S1 = sympy.Matrix([[sigma1,0],[0,sigma1]])
|
||||
S2 = sympy.Matrix([[sigma1,0],[0,sigma2]])
|
||||
S3 = sympy.Matrix([[sigma1,sigma3],[sigma3,sigma2]])
|
||||
|
||||
print '2D'
|
||||
print sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1))
|
||||
print sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1))
|
||||
print sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1))
|
||||
|
||||
y = 12
|
||||
z = 5
|
||||
ex = x**2+y*z
|
||||
e = ex
|
||||
|
||||
sigma1 = x*y+1
|
||||
|
||||
print '1D'
|
||||
print sympy.integrate(e*sigma1*e, (x,0,1))
|
||||
|
||||
Reference in New Issue
Block a user