From 435421731aa75b8ed2dc4020067308d30d54e23a Mon Sep 17 00:00:00 2001 From: carsten-forty2 Date: Wed, 20 Nov 2019 13:52:02 +0100 Subject: [PATCH] linting --- .pylintrc | 2 +- python/pygimli/solver/solver.py | 101 ++++++++++++++--------------- python/pygimli/testing/test_FEM.py | 38 +++++------ 3 files changed, 69 insertions(+), 72 deletions(-) diff --git a/.pylintrc b/.pylintrc index c65c5ee58..2326050fa 100644 --- a/.pylintrc +++ b/.pylintrc @@ -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). diff --git a/python/pygimli/solver/solver.py b/python/pygimli/solver/solver.py index e4ca52ab9..8663272ab 100644 --- a/python/pygimli/solver/solver.py +++ b/python/pygimli/solver/solver.py @@ -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 \ @@ -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. @@ -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 @@ -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) @@ -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 @@ -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) @@ -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. @@ -1202,13 +1204,13 @@ 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: @@ -1216,9 +1218,6 @@ def assembleBC_(bc, mesh, S, rhs, a, time=None, userData=None): 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. @@ -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:: @@ -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` @@ -1363,21 +1361,21 @@ 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 @@ -1385,7 +1383,7 @@ def H1Norm(u, S=None, mesh=None): u : iterable Node based value to compute the H1 norm for. - S : Matrix + mat : Matrix Stiffness matrix. mesh : :gimliapi:`GIMLI::Mesh` @@ -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): @@ -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') @@ -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) @@ -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: @@ -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)) @@ -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()) diff --git a/python/pygimli/testing/test_FEM.py b/python/pygimli/testing/test_FEM.py index 3f8343f01..047147dcc 100644 --- a/python/pygimli/testing/test_FEM.py +++ b/python/pygimli/testing/test_FEM.py @@ -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): """ @@ -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): """ """ @@ -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 @@ -106,19 +106,19 @@ 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 @@ -126,8 +126,8 @@ 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() @@ -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 @@ -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()]) @@ -182,7 +182,7 @@ def _testP2_(mesh, show=False): #TODO 3D Tet if __name__ == '__main__': - + # test = TestFiniteElementBasics() # test.test_Neumann() # test.test_Dirichlet()