Skip to content

Commit

Permalink
linting
Browse files Browse the repository at this point in the history
  • Loading branch information
@carsten-forty2
carsten-forty2 committed Nov 20, 2019
1 parent 10b0644 commit 4354217
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 72 deletions.
2 changes: 1 addition & 1 deletion .pylintrc
View file
Expand Up @@ -105,7 +105,7 @@ single-line-if-stmt=no
no-space-check=trailing-comma,dict-separator

# Maximum number of lines in a module
max-module-lines=1000
max-module-lines=2000

# String used as indentation unit. This is usually " " (4 spaces) or "\t" (1
# tab).
Expand Down
101 changes: 49 additions & 52 deletions python/pygimli/solver/solver.py
View file
Expand Up @@ -815,7 +815,8 @@ def linSolve(mat, b, solver=None, verbose=False):
if isinstance(mat, np.ndarray):
return np.linalg.solve(mat, b)

import scipy.sparse
scipy = pg.optImport('scipy')
#import scipy.sparse
from scipy.sparse.linalg import spsolve

if isinstance(mat, scipy.sparse.csr.csr_matrix) or \
Expand All @@ -829,17 +830,18 @@ def linSolve(mat, b, solver=None, verbose=False):

return x

def assembleForceVector(mesh, f, userData=None):
"""
DEPRECATED use assembleLoadVector instead
"""
return assembleLoadVector(mesh, f, userData)


def assembleLoadVector(mesh, f, userData=None):
r"""Assemble the load vector. See createLoadVector. Maybe we will remove this """
return createLoadVector(mesh, f, userData)

def createLoadVector(mesh, f, userData=None):
"""Create right hand side vector based on the given mesh and load
or force values.
TODO:
Maybe we can define an additional nodal force vector.
Create right hand side vector based on the given mesh and load or force
values.
Expand Down Expand Up @@ -869,7 +871,7 @@ def assembleLoadVector(mesh, f, userData=None):
rhs = np.zeros((len(f), mesh.nodeCount()))
for i, fi in enumerate(f):
userData['i'] = i
rhs[i] = assembleForceVector(mesh, fi, userData)
rhs[i] = createLoadVector(mesh, fi, userData)

return rhs

Expand Down Expand Up @@ -911,7 +913,7 @@ def assembleLoadVector(mesh, f, userData=None):
# for i, idx in enumerate(b_l.idx()):
# rhsRef[idx] += b_l.row(0)[i] * fArray[c.id()]
# np.testing.assert_allclose(rhs, rhsRef)
# print("Remove revtest in assembleForceVector after check")
# print("Remove revtest in assembleLoadVector after check")

elif len(fArray) == mesh.nodeCount():
fA = pg.Vector(fArray)
Expand All @@ -927,11 +929,11 @@ def assembleLoadVector(mesh, f, userData=None):
# for i, idx in enumerate(b_l.idx()):
# rhsRef[idx] += b_l.row(0)[i] * fArray[idx]
# np.testing.assert_allclose(rhs, rhsRef)
# print("Remove revtest in assembleForceVector after check")
# print("Remove revtest in assembleLoadVector after check")

# rhs = pg.Vector(fArray)
else:
raise Exception("Forcevector have the wrong size: " +
raise Exception("Load vector have the wrong size: " +
str(len(fArray)))

return rhs
Expand Down Expand Up @@ -1038,14 +1040,14 @@ def assembleDirichletBC(mat, boundaryPairs, rhs=None, time=0.0, userData=None,
if nodePairs is not None:
#print("nodePairs", nodePairs)

if len(nodePairs) == 2 and type(nodePairs[0]) == int:
if len(nodePairs) == 2 and isinstance(nodePairs[0], int):
# assume a single Node [NodeId, val]
nodePairs = [nodePairs]

for i, [n, val] in enumerate(nodePairs):
uDirIndex.append(n)
if hasattr(val, '__call__'):
raise("callabe node pairs need to be implement.")
raise "callabe node pairs need to be implement."
uDirichlet.append(val)

_assembleUDirichlet(mat, rhs, uDirIndex, uDirichlet)
Expand Down Expand Up @@ -1189,7 +1191,7 @@ def assembleRobinBC(mat, boundaryPairs, rhs=None, time=0.0, userData=None):
rhs.add(Sq, -p*q)


def assembleBC_(bc, mesh, S, rhs, a, time=None, userData=None):
def assembleBC_(bc, mesh, mat, rhs, a, time=None, userData=None):
r"""Shortcut to apply all boundary conditions.
This is a helper function for the solver call.
Expand All @@ -1202,23 +1204,20 @@ def assembleBC_(bc, mesh, S, rhs, a, time=None, userData=None):
assembleNeumannBC(rhs, parseArgToBoundaries(bct.pop('Neumann'), mesh),
a=a, time=time, userData=userData)
if 'Robin' in bct:
assembleRobinBC(S, parseArgToBoundaries(bct.pop('Robin'), mesh),
assembleRobinBC(mat, parseArgToBoundaries(bct.pop('Robin'), mesh),
rhs=rhs, time=time, userData=userData)
if 'Dirichlet' in bct:
assembleDirichletBC(S, parseArgToBoundaries(bct.pop('Dirichlet'), mesh),
assembleDirichletBC(mat, parseArgToBoundaries(bct.pop('Dirichlet'), mesh),
rhs=rhs, time=time, userData=userData)
if 'Node' in bct:
assembleDirichletBC(S, [], nodePairs=bct.pop('Node'),
assembleDirichletBC(mat, [], nodePairs=bct.pop('Node'),
rhs=rhs, time=time, userData=userData)

if len(bct.keys()) > 0:
pg.warn("Unknown boundary condition[s]" + \
str(bct.keys()) + " will be ignored")


def createLoadVector(mesh, f, userData=None):
return assembleLoadVector(mesh, f, userData)

def createStiffnessMatrix(mesh, a=None):
r"""Create the Stiffness matrix.
Expand Down Expand Up @@ -1318,19 +1317,18 @@ def createMassMatrix(mesh, b=None):
# return B


def _feNorm(u, A):
def _feNorm(u, mat):
"""Create a norm within a Finite Element space.
Create the Finite Element Norm with a preassembled system matrix.
"""
return np.sqrt(pg.math.dot(u, A.mult(u)))
return np.sqrt(pg.math.dot(u, mat.mult(u)))


def L2Norm(u, M=None, mesh=None):
r"""Create Lebesgue (L2) norm for the finite element space.
def normL2(u, mat=None, mesh=None):
r"""Create Lebesgue (L2) norm for finite element space.
Find the L2 Norm for a solution in the finite element space.
:math:`u` exact solution
Find the L2 Norm for a solution for the finite element space. :math:`u` exact solution
:math:`{\bf M}` Mass matrix, i.e., Finite element identity matrix.
.. math::
Expand All @@ -1351,7 +1349,7 @@ def L2Norm(u, M=None, mesh=None):
u : iterable
Node based value to compute the L2 norm for.
M : Matrix
mat : Matrix
Mass element matrix.
mesh : :gimliapi:`GIMLI::Mesh`
Expand All @@ -1363,29 +1361,29 @@ def L2Norm(u, M=None, mesh=None):
:math:`L2(u)` norm.
"""
if isinstance(M, pg.Mesh):
mesh = M
M = None
if isinstance(mat, pg.Mesh):
mesh = mat
mat = None

if M is None and mesh is not None:
M = createMassMatrix(mesh)
if mat is None and mesh is not None:
mat = createMassMatrix(mesh)

if M is None:
if mat is None:
# M is Identity matrix
return np.sqrt(pg.math.dot(u, u))

return _feNorm(u, M)
return _feNorm(u, mat)


def H1Norm(u, S=None, mesh=None):
def normH1(u, mat=None, mesh=None):
r"""Create (H1) norm for the finite element space.
Parameters
----------
u : iterable
Node based value to compute the H1 norm for.
S : Matrix
mat : Matrix
Stiffness matrix.
mesh : :gimliapi:`GIMLI::Mesh`
Expand All @@ -1397,17 +1395,17 @@ def H1Norm(u, S=None, mesh=None):
:math:`H1(u)` norm.
"""
if isinstance(S, pg.Mesh):
mesh = S
S = None
if isinstance(mat, pg.Mesh):
mesh = mat
mat = None

if S is None and mesh is not None:
S = pg.solver.createStiffnessMatrix(mesh)
if mat is None and mesh is not None:
mat = pg.solver.createStiffnessMatrix(mesh)

if S is None:
raise Exception("Need S or mesh here to calculate H1Norm")
if mat is None:
raise Exception("Need Stiffness matrix or a mesh here to calculate H1-Norm")

return _feNorm(u, S)
return _feNorm(u, mat)


def solve(mesh, **kwargs):
Expand Down Expand Up @@ -1534,7 +1532,7 @@ def solveFiniteElements(mesh, a=1.0, b=None, f=0.0, bc=None,
>>> plt.show()
"""
if bc is None:
bc={}
bc = {}
if 'uB' in kwargs:
pg.deprecated('bc arg uB', "bc={'Dirichlet': uB}") #20190912
bc['Dirichlet'] = kwargs.pop('uB')
Expand Down Expand Up @@ -1575,7 +1573,7 @@ def solveFiniteElements(mesh, a=1.0, b=None, f=0.0, bc=None,
A = S

if times is None:
rhs = assembleForceVector(mesh, f, userData=userData)
rhs = createLoadVector(mesh, f, userData=userData)

assembleBC_(bc, mesh, A, rhs, a, time=None, userData=userData)

Expand Down Expand Up @@ -1645,7 +1643,7 @@ def solveFiniteElements(mesh, a=1.0, b=None, f=0.0, bc=None,
print("start TL", swatch.duration())

M = createMassMatrix(mesh)
F = assembleForceVector(mesh, f)
F = assembleLoadVector(mesh, f)

u0 = np.zeros(dof)
if 'u0' in kwargs:
Expand Down Expand Up @@ -1724,10 +1722,9 @@ def solveFiniteElements(mesh, a=1.0, b=None, f=0.0, bc=None,
U[n, :] = np.asarray(u)

if progress:
progress.update(n,
' t_prep: ' + str(round(t_prep, 5)) + 's' + \
' t_step: ' + str(round(swatch.duration(),5)) + 's')

progress.update(n, ' t_prep: {0}ms t_step {1}s'.format(
pg.pf(t_prep*1000),
pg.pf(swatch.duration())))
if debug:
print("Measure(" + str(len(times)) + "): ",
measure, measure / len(times))
Expand Down Expand Up @@ -1828,7 +1825,7 @@ def crankNicolson(times, theta, S, I, f, u0=None, progress=None):
timeAssemble.append(pg.dur())
pg.tic()

u[n,:] = solver.solve(b)
u[n, :] = solver.solve(b)

if timeMeasure:
timeSolve.append(pg.dur())
Expand Down
38 changes: 19 additions & 19 deletions python/pygimli/testing/test_FEM.py
View file
Expand Up @@ -19,13 +19,13 @@ def test_Poisson(self):
f = -cos(x)
### the following are already tested in Helmholtz
u = x
u = x
f = 0
u = x*x
f = -2
"""
pass
pass

def test_Helmholtz(self):
"""
Expand Down Expand Up @@ -65,7 +65,7 @@ def test_Helmholtz(self):
# pg.plt.plot(x, uFEM, '.')
# pg.plt.plot(pg.sort(x), u(pg.sort(x)))
# pg.wait()

def test_Neumann(self):
"""
"""
Expand All @@ -82,7 +82,7 @@ def _test_(mesh, show=False):
pg.show(grid, pg.abs(v))
pg.show(grid, v, showMesh=1)
pg.wait()

v = pg.solver.grad(mesh, u)
np.testing.assert_allclose(pg.abs(v), np.ones(mesh.cellCount())*vTest)
return v
Expand All @@ -106,28 +106,28 @@ def _testP1_(mesh, show=False):
""" Laplace u = 0 solves u = x for u(r=0)=0 and u(r=1)=1
Test for u == exact x for P1 base functions
"""
u = pg.solve(mesh, a=1, b=0, f=0,
u = pg.solve(mesh, a=1, b=0, f=0,
bc={'Dirichlet': [[1, 0], [2, 1]]})

if show:
if mesh.dim()==1:
if mesh.dim()==1:
pg.plt.plot(pg.x(mesh), u)
pg.wait()
elif mesh.dim()==2:
pg.show(mesh, u, label='u')
pg.wait()

xMin = mesh.xmin()
xSpan = (mesh.xmax() - xMin)
xMin = mesh.xMin()
xSpan = (mesh.xMax() - xMin)
np.testing.assert_allclose(u, (pg.x(mesh)-xMin) / xSpan)
return u

def _testP2_(mesh, show=False):
""" Laplace u = 2 solves u = x² for u(r=0)=0 and u(r=1)=1
Test for u == exact x² for P2 base functions
"""
meshp2 = mesh.createP2()
u = pg.solve(meshp2, f=-2, bc={'Dirichlet': [[1, 0], [2, 1]]})
meshP2 = mesh.createP2()
u = pg.solve(meshP2, f=-2, bc={'Dirichlet': [[1, 0], [2, 1]]})

# find test pos different from node pos
meshTests = mesh.createH2()
Expand All @@ -140,22 +140,22 @@ def _testP2_(mesh, show=False):
c = [b.center() for b in meshTests.boundaries(meshTests.boundaryMarkers()==4)]

c.sort(key=lambda c_: c_.distance(startPos))
ui = pg.interpolate(meshp2, u, c)
xi = pg.utils.cumDist(c) + startPos.distance(c[0])
ui = pg.interpolate(meshP2, u, c)
xi = pg.utils.cumDist(c) + startPos.distance(c[0])

if show:
pg.plt.plot(xi, ui)
pg.plt.plot(xi, xi**2)
pg.wait()

np.testing.assert_allclose(ui, xi**2)

_testP1_(pg.createGrid(x=np.linspace(0, 1, 11)), show=False) #1D
_testP1_(pg.createGrid(x=np.linspace(-2, 1, 11),
_testP1_(pg.createGrid(x=np.linspace(-2, 1, 11),
y=np.linspace(0, 1, 11))) #2D reg quad
_testP1_(pg.createGrid(x=np.linspace(-0.04, 0.01, 11),
_testP1_(pg.createGrid(x=np.linspace(-0.04, 0.01, 11),
y=np.linspace(-0.4, 0, 11))) #2D scaled
_testP1_(pg.createGrid(x=np.linspace(-2, 1, 11),
_testP1_(pg.createGrid(x=np.linspace(-2, 1, 11),
y=np.linspace( 0, 1, 11),
z=np.linspace( 0, 1, 11))) #3D

Expand All @@ -167,7 +167,7 @@ def _testP2_(mesh, show=False):
grid.rotate([0, 0, np.pi/4])
#can't find proper test for this rotated case
#v = _testP1_(grid, show=False) #2D reg - rotated

# P2 tests
mesh = pg.meshtools.createMesh(pg.meshtools.createWorld(start=[0, -4], end=[1, -6], worldMarker=0), area=0.1)
mesh.setBoundaryMarkers(np.array([0,1,3,2,4])[mesh.boundaryMarkers()])
Expand All @@ -182,7 +182,7 @@ def _testP2_(mesh, show=False):
#TODO 3D Tet

if __name__ == '__main__':

# test = TestFiniteElementBasics()
# test.test_Neumann()
# test.test_Dirichlet()
Expand Down

0 comments on commit 4354217

Please sign in to comment.