parameters: (mesh, a=1.0, b=None, f=0.0, bc=None,
                          times=None, c=1.0, userData={},
                          verbose=False, **kwargs)
docstring: r"""Solve partial differential equation with Finite Elements.
  
      This is a syntactic sugar convenience function for using the Finite Element
      functionality of the library core to solve partial differential equation
      (PDE) that match the following form:
  
      .. math::
  
          c \frac{\partial u}{\partial t} & = \nabla\cdot(a \nabla u)
          + b u + f(\mathbf{r},t)~~|~~\Omega_{\text{Mesh}}\\
          u & = h~~|~~\Gamma_{\text{Dirichlet}}\\
          \frac{\partial u}{\partial \mathbf{n}} & = g~~|~~\Gamma_{\text{Neumann}}\\
          \alpha u + \beta\frac{\partial u}{\partial \mathbf{n}} & = \gamma~~|~~\Gamma_{\text{Robin}}\\
          \frac{\partial u}{\partial \mathbf{n}} & = \alpha(u_0-u)~~|~~\Gamma_{\text{Robin}}
  
      for the scalar :math:`u(\mathbf{r}, t)` or
      vector :math:`\mathbf(u)(\mathbf{r}, t)` solution at each node of a given mesh.
      The Domain :math:`\Omega` and the Boundary :math:`\Gamma` are defined
      through the mesh with appropriate boundary marker.
  
      Note, to ensure vector solution either set vector forces or at least on
      vector component boundary condition.
  
      TODO:
  
          * unsteady ub and dub
          * 'Infinity' Boundary condition (u vanishes at infinity)
          * 'Cauchy' Boundary condition
          (guaranties u and du on same boundary, will never work here because the
          problem becomes ill posed and would need some inverse strategy to
          solve.)
          * Example for
              * elastic parameter
              * anisotropic (float/complex)
              * dynamic boundary conditions
              * dynamic load vector
              * nonlinearity
  
      Parameters
      ----------
      mesh: :gimliapi:`GIMLI::Mesh`
          Mesh represents spatial discretization of the calculation domain
      a: value | array | callable(cell, userData)
          Cell values of type float or complex can be scalar, anisotropic matrix or elastic tensor.
      b: value | array | callable(cell, userData) [None]
          Cell values. None means the term is unused.
      c: value | array | callable(cell, userData) [None]
          Scale the unsteady term, only for times is not None.
      f: value | array(cells) | array(nodes) | callable(args, kwargs)
          force values, for vector fields use (n x dim) values.
      bc: dict()
          Dictionary of boundary conditions.
          Current supported boundary conditions by dictionary keys:
          'Dirichlet', 'Neumann', 'Robin', 'Node'.
  
          The dictionary can contain multiple "key: Arg"
          Arg will be parsed by :py:mod:`pygimli.solver.solver.parseArgPairToBoundaryArray`
  
          If the dictionary key is 'Node' then fixed values for single node
          indices can by be given. e.g., bc={'Node': [nodeID, value]}.
          Note this is only a shortcut for
          bc={'Dirichlet': [mesh.node(nodeID), value]}.
      times: array [None]
          Solve as time dependent problem for the given times.
  
      Keyword Arguments
      -----------------
      **kwargs
          u0: value | array | callable(pos, userData)
              Node values
          theta: float [1]
              - :math:`theta = 0` means explicit Euler, maybe stable for
              :math:`\Delta t \quad\text{near}\quad h`
              - :math:`theta = 0.5`, Crank-Nicolson scheme, maybe instable
              - :math:`theta = 2/3`, Galerkin scheme
              - :math:`theta = 1`, implicit Euler
  
              If unsure choose :math:`\theta = 0.5 + \epsilon`, which is probably stable.
          dynamic: bool [False]
              Boundary conditions for time depending problems will be considered
              dynamic for each time step.
          stats: bool
              Give some statistics.
          progress: bool
              Give some calculation progress.
          assembleOnly: bool
              Stops after matrix asssemblation.
              Returns the system matrix A and the rhs vector.
          fixPureNeumann: bool [auto]
              If set or detected automatic, we add the additional condition:
              :math:`\int_domain u dv = 0` which makes elliptic problems well posed again.
          rhs: iterable
              Pre assembled rhs. Will preferred on any f settings.
          ws: dict
              The WorkSpace is a dictionary that will get
              some temporary data during the calculation.
              Any keyvalue 'u' in the dictionary will be used for the resulting array.
          vectorValued: bool (False)
              Solution forced to vector valued, in case the auto detection fails
  
      Returns
      -------
      u: array
          Returns the solution u either 1,n array for stationary problems or
          for m,n array for m time steps
  
      See also
      --------
          :ref:`tut:modelling` and :py:mod:`pygimli.solver.solve`
  
      Examples
      --------
      >>> import pygimli as pg
      >>> from pygimli.meshtools import polytools as plc
      >>> from pygimli.viewer.mpl import drawField, drawMesh
      >>> import matplotlib.pyplot as plt
      >>> world = plc.createWorld(start=[-10, 0], end=[10, -10],
      ...                         marker=1, worldMarker=False)
      >>> c1 = plc.createCircle(pos=[0.0, -5.0], radius=3.0, area=.1, marker=2)
      >>> mesh = pg.meshtools.createMesh([world, c1], quality=34.3)
      >>> u = pg.solver.solveFiniteElements(mesh, a={1: 100.0, 2: 1.0},
      ...                                   bc={'Dirichlet':{4: 1.0, 2: 0.0}})
      >>> fig, ax = plt.subplots()
      >>> pc = drawField(ax, mesh, u)
      >>> drawMesh(ax, mesh)
      >>> plt.show()
docstring_summary: r"""Solve partial differential equation with Finite Elements.
docstring_tokens: r Solve partial differential equation with Finite Elements .
function: def solveFiniteElements(mesh, a=1.0, b=None, f=0.0, bc=None,
                          times=None, c=1.0, userData={},
                          verbose=False, **kwargs):
      r"""Solve partial differential equation with Finite Elements.
  
      This is a syntactic sugar convenience function for using the Finite Element
      functionality of the library core to solve partial differential equation
      (PDE) that match the following form:
  
      .. math::
  
          c \frac{\partial u}{\partial t} & = \nabla\cdot(a \nabla u)
          + b u + f(\mathbf{r},t)~~|~~\Omega_{\text{Mesh}}\\
          u & = h~~|~~\Gamma_{\text{Dirichlet}}\\
          \frac{\partial u}{\partial \mathbf{n}} & = g~~|~~\Gamma_{\text{Neumann}}\\
          \alpha u + \beta\frac{\partial u}{\partial \mathbf{n}} & = \gamma~~|~~\Gamma_{\text{Robin}}\\
          \frac{\partial u}{\partial \mathbf{n}} & = \alpha(u_0-u)~~|~~\Gamma_{\text{Robin}}
  
      for the scalar :math:`u(\mathbf{r}, t)` or
      vector :math:`\mathbf(u)(\mathbf{r}, t)` solution at each node of a given mesh.
      The Domain :math:`\Omega` and the Boundary :math:`\Gamma` are defined
      through the mesh with appropriate boundary marker.
  
      Note, to ensure vector solution either set vector forces or at least on
      vector component boundary condition.
  
      TODO:
  
          * unsteady ub and dub
          * 'Infinity' Boundary condition (u vanishes at infinity)
          * 'Cauchy' Boundary condition
          (guaranties u and du on same boundary, will never work here because the
          problem becomes ill posed and would need some inverse strategy to
          solve.)
          * Example for
              * elastic parameter
              * anisotropic (float/complex)
              * dynamic boundary conditions
              * dynamic load vector
              * nonlinearity
  
      Parameters
      ----------
      mesh: :gimliapi:`GIMLI::Mesh`
          Mesh represents spatial discretization of the calculation domain
      a: value | array | callable(cell, userData)
          Cell values of type float or complex can be scalar, anisotropic matrix or elastic tensor.
      b: value | array | callable(cell, userData) [None]
          Cell values. None means the term is unused.
      c: value | array | callable(cell, userData) [None]
          Scale the unsteady term, only for times is not None.
      f: value | array(cells) | array(nodes) | callable(args, kwargs)
          force values, for vector fields use (n x dim) values.
      bc: dict()
          Dictionary of boundary conditions.
          Current supported boundary conditions by dictionary keys:
          'Dirichlet', 'Neumann', 'Robin', 'Node'.
  
          The dictionary can contain multiple "key: Arg"
          Arg will be parsed by :py:mod:`pygimli.solver.solver.parseArgPairToBoundaryArray`
  
          If the dictionary key is 'Node' then fixed values for single node
          indices can by be given. e.g., bc={'Node': [nodeID, value]}.
          Note this is only a shortcut for
          bc={'Dirichlet': [mesh.node(nodeID), value]}.
      times: array [None]
          Solve as time dependent problem for the given times.
  
      Keyword Arguments
      -----------------
      **kwargs
          u0: value | array | callable(pos, userData)
              Node values
          theta: float [1]
              - :math:`theta = 0` means explicit Euler, maybe stable for
              :math:`\Delta t \quad\text{near}\quad h`
              - :math:`theta = 0.5`, Crank-Nicolson scheme, maybe instable
              - :math:`theta = 2/3`, Galerkin scheme
              - :math:`theta = 1`, implicit Euler
  
              If unsure choose :math:`\theta = 0.5 + \epsilon`, which is probably stable.
          dynamic: bool [False]
              Boundary conditions for time depending problems will be considered
              dynamic for each time step.
          stats: bool
              Give some statistics.
          progress: bool
              Give some calculation progress.
          assembleOnly: bool
              Stops after matrix asssemblation.
              Returns the system matrix A and the rhs vector.
          fixPureNeumann: bool [auto]
              If set or detected automatic, we add the additional condition:
              :math:`\int_domain u dv = 0` which makes elliptic problems well posed again.
          rhs: iterable
              Pre assembled rhs. Will preferred on any f settings.
          ws: dict
              The WorkSpace is a dictionary that will get
              some temporary data during the calculation.
              Any keyvalue 'u' in the dictionary will be used for the resulting array.
          vectorValued: bool (False)
              Solution forced to vector valued, in case the auto detection fails
  
      Returns
      -------
      u: array
          Returns the solution u either 1,n array for stationary problems or
          for m,n array for m time steps
  
      See also
      --------
          :ref:`tut:modelling` and :py:mod:`pygimli.solver.solve`
  
      Examples
      --------
      >>> import pygimli as pg
      >>> from pygimli.meshtools import polytools as plc
      >>> from pygimli.viewer.mpl import drawField, drawMesh
      >>> import matplotlib.pyplot as plt
      >>> world = plc.createWorld(start=[-10, 0], end=[10, -10],
      ...                         marker=1, worldMarker=False)
      >>> c1 = plc.createCircle(pos=[0.0, -5.0], radius=3.0, area=.1, marker=2)
      >>> mesh = pg.meshtools.createMesh([world, c1], quality=34.3)
      >>> u = pg.solver.solveFiniteElements(mesh, a={1: 100.0, 2: 1.0},
      ...                                   bc={'Dirichlet':{4: 1.0, 2: 0.0}})
      >>> fig, ax = plt.subplots()
      >>> pc = drawField(ax, mesh, u)
      >>> drawMesh(ax, mesh)
      >>> plt.show()
      """
      if bc is None:
          bc = {}
  
      workSpace = kwargs.pop('ws', dict())
      debug = kwargs.pop('debug', False)
      stats = kwargs.pop('stats', False)
  
      mesh.createNeighborInfos()
      if verbose:
          print("Mesh: ", str(mesh))
  
      ## scalar solution default
      vectorValues = False
      dof = mesh.nodeCount()
  
      ## check if force vector is a vector
      rhs = kwargs.pop('rhs', createLoadVector(mesh, f, userData=userData))
  
      # pg._g('###############')
      if len(rhs) > dof or kwargs.pop('vectorValued',
                                      _bcIsForVectorValues(bc, mesh)):
          if verbose:
              print("Solve vector valued.")
          vectorValues = True
          dof = mesh.nodeCount() * mesh.dimension()
  
      # pg._g('###############', dof)
      rhs.resize(dof)
  
      swatch = pg.core.Stopwatch(True)
  
      # check for material parameter
      #a = parseArgToArray(a, nDof=mesh.cellCount(), mesh=mesh, userData=userData)
      a = cellValues(mesh, a, userData=userData)
      isComplex = False
      if pg.utils.isComplex(a):
          isComplex = True
          rhs = np.array(rhs, dtype=complex)
  
      S = createStiffnessMatrix(mesh, a, isVector=vectorValues)
      M = None
  
      if b is not None and b != 0:
          b = cellValues(mesh, b, userData=userData)
          M = createMassMatrix(mesh, b)
          #pg.warn("check me")
          A = S - M
      else:
          A = S
  
      if times is None:
  
          if len(list(bc.items())) == 0 or \
             (len(list(bc.items())) == 1 and list(bc.keys())[0] == 'Neumann'):
              pn = True
          else:
              pn = False
  
          fixPureNeumann = kwargs.pop('fixPureNeumann', pn)
  
          assembleBC(bc, mesh, A, rhs, time=None, userData=userData)
  
          u = None
          if 'u' in workSpace:
              u = workSpace['u']
  
          singleForce = True
  
          if hasattr(rhs, 'ndim'):
              if rhs.ndim == 2:
                  singleForce = False
                  if u is None:
                      u = np.zeros(rhs.shape)
          else:
              if isinstance(a[0], complex):
                  if u is None:
                      u = pg.CVector(rhs.size(), 0.0)
                  rhs = pg.math.toComplex(rhs)
              else:
                  if u is None:
                      u = pg.Vector(rhs.size(), 0.0)
  
          if fixPureNeumann is True:
              pg.info('Fixing pure Neumann boundary condition by forcing: '
                      'intDomain(u, mesh) = 0')
              r = createLoadVector(mesh)
  
              A = pg.BlockMatrix()
              A.add(S, 0, 0)
              A.add(r, 0, mesh.nodeCount())
              A.add(r, mesh.nodeCount(), 0, transpose=True)
  
              rhs = pg.cat(rhs, pg.Vector(1,0))
  
          assembleTime = swatch.duration(True)
  
          if stats:
              stats.assembleTime = assembleTime
  
          if verbose:
              print("Assembling time: ", assembleTime)
  
          workSpace['Stiffness matrix'] = S
          workSpace['Mass matrix'] = M
          workSpace['System matrix'] = A
          workSpace['rhs'] = rhs
  
          if 'assembleOnly' in kwargs:
              return A, rhs
  
          if fixPureNeumann:
              if singleForce:
                  uc = pg.solver.linSolve(A, rhs, 'scipy')
                  u = uc[0:mesh.nodeCount()]
              else:
                  pg.criticale('Non single force for pure Neumann not yet implemented')
          else:
              solver = pg.core.LinSolver(False)
              solver.setMatrix(A, 0)
  
              if singleForce:
                  if isComplex is True:
                      # clean this up
                      rhs = pg.core.toComplex(rhs.real, rhs.imag)
                      u = solver.solve(rhs).array()
                  else:
                      u = solver.solve(rhs)
              else:
                  for i, r in enumerate(rhs):
                      u[i] = solver.solve(r)
  
          solverTime = swatch.duration(True)
          if verbose:
              if stats:
                  stats.solverTime = solverTime
              print("Solving time: ", solverTime)
  
          if len(kwargs.keys()) > 0:
              pg.warn("Unused arguments", *kwargs)
          return u
  
      else: # times given
          pg.solver.checkCFL(times, mesh, max(a))
  
          if debug:
              print("start TL", swatch.duration())
  
          if c != 1.0:
              c = cellValues(mesh, c, userData=userData)
  
          M = createMassMatrix(mesh, c)
          F = createLoadVector(mesh, f)
  
          u0 = np.zeros(dof)
          if 'u0' in kwargs:
              u0 = parseArgToArray(kwargs['u0'], dof, mesh, userData)
  
          progress = None
          if 'progress' in kwargs:
              from pygimli.utils import ProgressBar
              progress = ProgressBar(its=len(times), width=40, sign='+')
  
          theta = kwargs.pop('theta', 1.0)
          dynamic = kwargs.pop('dynamic', False)
  
          if not dynamic:
              S = createStiffnessMatrix(mesh, a)
              assembleBC(bc, mesh, S, F, time=0.0, userData=userData)
              return crankNicolson(times, S, M, f=F,
                                   u0=u0, theta=theta,
                                   progress=progress)
  
          rhs = np.zeros((len(times), dof))
          # no time dependency for rhs so far ... TODO
          rhs[:] = F  # this is slow: optimize
  
          if debug:
              print("rhs", swatch.duration())
          U = np.zeros((len(times), dof))
          U[0, :] = u0
  
          # init state
          u = pg.Vector(dof, 0.0)
  
          if debug:
              print("u0", swatch.duration())
  
          measure = 0.
          for n in range(1, len(times)):
              swatch.reset()
  
              dt = times[n] - times[n-1]
  
              # previous timestep
              # print "i: ", i, dt, U[i - 1]
  
              swatch.reset()
              # (A + a*B)u is fastest,
              # followed by A*u + (B*u)*a and finally A*u + a*B*u and
              br = (M + (dt * (theta - 1.)) * S) * U[n - 1] + \
                  dt * ((1.0 - theta) * rhs[n - 1] + theta * rhs[n])
  
              # print ('a',swatch.duration(True))
              # br = M * U[n - 1] - (A * U[n - 1]) * (dt*(1.0 - theta)) + \
              # dt * ((1.0 - theta) * rhs[n - 1] + theta * rhs[n])
  
              # print ('br',swatch.duration(True))
  
              # br = M * U[n - 1] - (dt*(1.0 - theta)) * A * U[n - 1] + \
              # dt * ((1.0 - theta) * rhs[n - 1] + theta * rhs[n])
              # print ('c',swatch.duration(True))
  
              measure += swatch.duration()
  
              A = M + S * dt * theta
  
              assembleBC(bc, mesh, A, br, time=times[n], userData=userData)
  
              if 'assembleOnly' in kwargs:
                  return A, br
  
              # u = S/b
              t_prep = swatch.duration(True)
              solver = pg.core.LinSolver(A, verbose)
              solver.solve(br, u)
  
              if 'plotTimeStep' in kwargs:
                  kwargs['plotTimeStep'](u, times[n])
  
              U[n, :] = np.asarray(u)
  
              if progress:
                  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))
          return U
function_tokens: def solveFiniteElements ( mesh , a = 1.0 , b = None , f = 0.0 , bc = None , times = None , c = 1.0 , userData = { } , verbose = False , * * kwargs ) : if bc is None : bc = { } workSpace = kwargs . pop ( 'ws' , dict ( ) ) debug = kwargs . pop ( 'debug' , False ) stats = kwargs . pop ( 'stats' , False ) mesh . createNeighborInfos ( ) if verbose : print ( "Mesh: " , str ( mesh ) ) ## scalar solution default vectorValues = False dof = mesh . nodeCount ( ) ## check if force vector is a vector rhs = kwargs . pop ( 'rhs' , createLoadVector ( mesh , f , userData = userData ) ) # pg._g('###############') if len ( rhs ) > dof or kwargs . pop ( 'vectorValued' , _bcIsForVectorValues ( bc , mesh ) ) : if verbose : print ( "Solve vector valued." ) vectorValues = True dof = mesh . nodeCount ( ) * mesh . dimension ( ) # pg._g('###############', dof) rhs . resize ( dof ) swatch = pg . core . Stopwatch ( True ) # check for material parameter #a = parseArgToArray(a, nDof=mesh.cellCount(), mesh=mesh, userData=userData) a = cellValues ( mesh , a , userData = userData ) isComplex = False if pg . utils . isComplex ( a ) : isComplex = True rhs = np . array ( rhs , dtype = complex ) S = createStiffnessMatrix ( mesh , a , isVector = vectorValues ) M = None if b is not None and b != 0 : b = cellValues ( mesh , b , userData = userData ) M = createMassMatrix ( mesh , b ) #pg.warn("check me") A = S - M else : A = S if times is None : if len ( list ( bc . items ( ) ) ) == 0 or ( len ( list ( bc . items ( ) ) ) == 1 and list ( bc . keys ( ) ) [ 0 ] == 'Neumann' ) : pn = True else : pn = False fixPureNeumann = kwargs . pop ( 'fixPureNeumann' , pn ) assembleBC ( bc , mesh , A , rhs , time = None , userData = userData ) u = None if 'u' in workSpace : u = workSpace [ 'u' ] singleForce = True if hasattr ( rhs , 'ndim' ) : if rhs . ndim == 2 : singleForce = False if u is None : u = np . zeros ( rhs . shape ) else : if isinstance ( a [ 0 ] , complex ) : if u is None : u = pg . CVector ( rhs . size ( ) , 0.0 ) rhs = pg . math . toComplex ( rhs ) else : if u is None : u = pg . Vector ( rhs . size ( ) , 0.0 ) if fixPureNeumann is True : pg . info ( 'Fixing pure Neumann boundary condition by forcing: ' 'intDomain(u, mesh) = 0' ) r = createLoadVector ( mesh ) A = pg . BlockMatrix ( ) A . add ( S , 0 , 0 ) A . add ( r , 0 , mesh . nodeCount ( ) ) A . add ( r , mesh . nodeCount ( ) , 0 , transpose = True ) rhs = pg . cat ( rhs , pg . Vector ( 1 , 0 ) ) assembleTime = swatch . duration ( True ) if stats : stats . assembleTime = assembleTime if verbose : print ( "Assembling time: " , assembleTime ) workSpace [ 'Stiffness matrix' ] = S workSpace [ 'Mass matrix' ] = M workSpace [ 'System matrix' ] = A workSpace [ 'rhs' ] = rhs if 'assembleOnly' in kwargs : return A , rhs if fixPureNeumann : if singleForce : uc = pg . solver . linSolve ( A , rhs , 'scipy' ) u = uc [ 0 : mesh . nodeCount ( ) ] else : pg . criticale ( 'Non single force for pure Neumann not yet implemented' ) else : solver = pg . core . LinSolver ( False ) solver . setMatrix ( A , 0 ) if singleForce : if isComplex is True : # clean this up rhs = pg . core . toComplex ( rhs . real , rhs . imag ) u = solver . solve ( rhs ) . array ( ) else : u = solver . solve ( rhs ) else : for i , r in enumerate ( rhs ) : u [ i ] = solver . solve ( r ) solverTime = swatch . duration ( True ) if verbose : if stats : stats . solverTime = solverTime print ( "Solving time: " , solverTime ) if len ( kwargs . keys ( ) ) > 0 : pg . warn ( "Unused arguments" , * kwargs ) return u else : # times given pg . solver . checkCFL ( times , mesh , max ( a ) ) if debug : print ( "start TL" , swatch . duration ( ) ) if c != 1.0 : c = cellValues ( mesh , c , userData = userData ) M = createMassMatrix ( mesh , c ) F = createLoadVector ( mesh , f ) u0 = np . zeros ( dof ) if 'u0' in kwargs : u0 = parseArgToArray ( kwargs [ 'u0' ] , dof , mesh , userData ) progress = None if 'progress' in kwargs : from pygimli . utils import ProgressBar progress = ProgressBar ( its = len ( times ) , width = 40 , sign = '+' ) theta = kwargs . pop ( 'theta' , 1.0 ) dynamic = kwargs . pop ( 'dynamic' , False ) if not dynamic : S = createStiffnessMatrix ( mesh , a ) assembleBC ( bc , mesh , S , F , time = 0.0 , userData = userData ) return crankNicolson ( times , S , M , f = F , u0 = u0 , theta = theta , progress = progress ) rhs = np . zeros ( ( len ( times ) , dof ) ) # no time dependency for rhs so far ... TODO rhs [ : ] = F # this is slow: optimize if debug : print ( "rhs" , swatch . duration ( ) ) U = np . zeros ( ( len ( times ) , dof ) ) U [ 0 , : ] = u0 # init state u = pg . Vector ( dof , 0.0 ) if debug : print ( "u0" , swatch . duration ( ) ) measure = 0. for n in range ( 1 , len ( times ) ) : swatch . reset ( ) dt = times [ n ] - times [ n - 1 ] # previous timestep # print "i: ", i, dt, U[i - 1] swatch . reset ( ) # (A + a*B)u is fastest, # followed by A*u + (B*u)*a and finally A*u + a*B*u and br = ( M + ( dt * ( theta - 1. ) ) * S ) * U [ n - 1 ] + dt * ( ( 1.0 - theta ) * rhs [ n - 1 ] + theta * rhs [ n ] ) # print ('a',swatch.duration(True)) # br = M * U[n - 1] - (A * U[n - 1]) * (dt*(1.0 - theta)) + \ # dt * ((1.0 - theta) * rhs[n - 1] + theta * rhs[n]) # print ('br',swatch.duration(True)) # br = M * U[n - 1] - (dt*(1.0 - theta)) * A * U[n - 1] + \ # dt * ((1.0 - theta) * rhs[n - 1] + theta * rhs[n]) # print ('c',swatch.duration(True)) measure += swatch . duration ( ) A = M + S * dt * theta assembleBC ( bc , mesh , A , br , time = times [ n ] , userData = userData ) if 'assembleOnly' in kwargs : return A , br # u = S/b t_prep = swatch . duration ( True ) solver = pg . core . LinSolver ( A , verbose ) solver . solve ( br , u ) if 'plotTimeStep' in kwargs : kwargs [ 'plotTimeStep' ] ( u , times [ n ] ) U [ n , : ] = np . asarray ( u ) if progress : 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 ) ) return U