evalf/nutils

The nutils project

nutils's Issues

Gmesh reader physical vs elementary groups

The gmesh reader is not adequately dealing with the group IDs in the case that physical entities are employed. From the gmesh manual:

number-of-tags
gives the number of integer tags that follow for the n-th element. By default, the first tag is the number of the physical entity to which the element belongs; the second is the number of the elementary geometrical entity to which the element belongs; the third is the number of mesh partitions to which the element belongs, followed by the partition ids (negative partition ids indicate ghost cells). A zero tag is equivalent to no tag. Gmsh and most codes using the MSH 2 format require at least the first two tags (physical and elementary tags).

Non-integer problem in rational.py

A machine dependent error occurs in the rational module where a numpy.int64 is returned as the numerator, whereas a int is expected.

Unittests fail for numpy version < 1.8(.1)

Many of the unittests fail when using numpy versions prior to 1.8. The suggested "solution" is to check for the numpy version in the install script and force the user to upgrade if his/her version is too old.

prepare for neumann

In preparing a spline basis for enforcing Neumann conditions strongly, one has at least two options:
a) h0 = h0+h1
b) h1 = h0+h1
In previous versions (finity) option (b) was chosen.

I would like to make a case for option (a). In both cases a Neumann constraint would affect only one DOF and avoid tedious linear constraints. But only option (a) recovers these advantages in the case that one would like to impose a Dirichlet constraint there too. I foresee this happening with higher order PDEs, or for the sake of recycling, or for care-free toggling between different conditions.

progress() and nprocs>1

Is it an idea to have progress() print only if pid=0? Printing progress for all procs does not give more info, right?

problem with @fork on some machines

Today I noticed a problem with the @fork decorator which causes the function to hang on my atom processor laptop; for instance the plotting routine in the laplace example script. The problem does not occur on mefd02, and it did not use to occur on my laptop, either, though I cannot determine what changed in the meantime.

I managed to trace the problem to numpy's lapack wrappers, which are inside numpy.linalg.lapack_lite. The following code should print 'done' twice, but on my laptop does so only once. The forked process hangs at 100% cpu.

import os, numpy

def main():
  pid = os.fork()
  if pid:
    print 'forked', pid
  else:
    print numpy.linalg.inv( numpy.eye( 2 ) )
  print 'done'

main()

Let me know if this bug affects you as well. Removing the fork decorator serves as a temporary workaround until I figure out a proper solution.

add support for complex arrays, functions

This enables solving the Helmholtz equation.

Signed up

I am watching the project on GitHub now.

StructuredTopology.linearfunc

Does line 955 in topology.py of the current commit (72c7d51) purposely disable an otherwise perfectly usable linearfunc() inherited from its parent class?

DeprecationWarning: `scipy.sparse.sparsetools` is deprecated!

If you get the following warning:

DeprecationWarning: scipy.sparse.sparsetools is deprecated!
scipy.sparse.sparsetools is a private module for scipy.sparse, and should not be used.

then maybe you're working with a newer version of scipy (like me, I'm using 0.14.0). A dirty workaround is changing the import in matrix.py to

from scipy.sparse import csr as _csr

But this is dirty. There must be a better way? Probably Trilinos will save us eventually from using scipy at all?

random errors drivencavity unittest

As of commit 5a3db07 the drivencavity unittest fails 1 out of 2 times when using python3.4 on debian stretch amd64 and fails always using python2.7 and python3.4 on debian stretch armhf.

On armhf testing commit c207d13:

$ python examples/drivencavity unittest --tbexplore=False --verbose=3
checkdata > compare > #1 > non matching lenghts: 6 != 29
checkdata > compare > #1 > #2 > non equal to 4 digits: 4.369273e-14 != 2.900000e-14
checkdata > compare > #1 > #3 > non equal to 4 digits: 3.309304e-14 != 1.000000e-15
checkdata > compare > #1 > #4 > non equal to 4 digits: 3.305709e-14 != 1.000000e-15
checkdata > compare > #1 > #5 > non equal to 4 digits: 2.115575e-14 != 0.000000e+00
checkdata > objects are not equal; if this is expected replace base64 data with:
eNqNkl1uAzEIhK+TSHbF8GsOlPtfoTbstnlJ1RezgvH4Y7QYDx2w53g8pvGy17ShQXLqhCS/pg5XidMA
HYGOyUlUCoZod9ai88FBOFUTqwbKmi3NbCk8S4q0ktiKcnWgTW0BfWW5n4mYF5BZlPkNOhVcSpYo0B9A
4JDvgXZDMrTfRDl66L7oI9kvOu8VNVY0BOLmFrRn2LWAoUa/wIYOzHm15CaFd1C7lgKhWlzk1XdK77lF
97m4Z4TwAZw7vnpMXFY7QPtGO+0MoolTljcFqFLZy+YllStP0jKVK56ZeVbwLcl2hTZPdrDCQh2wVcD0
NfavgpDygZwKjIOCrZb3chw+FbQE+KPQe9EPp/zj5LezvZ/j+Q14zqVr

AssertionError()
  File "examples/drivencavity", line 84, in unittest
    qe0DvRUYct7B2fr328qvHwXKg10=''' )
  File "nutils/util.py", line 524, in run
    func( **kwargs )

$ python3 examples/drivencavity unittest --tbexplore=False --verbose=3
checkdata > compare > #1 > non matching lenghts: 6 != 29
checkdata > compare > #1 > #2 > non equal to 4 digits: 4.369273e-14 != 2.900000e-14
checkdata > compare > #1 > #3 > non equal to 4 digits: 3.309304e-14 != 1.000000e-15
checkdata > compare > #1 > #4 > non equal to 4 digits: 3.305709e-14 != 1.000000e-15
checkdata > compare > #1 > #5 > non equal to 4 digits: 2.115575e-14 != 0.000000e+00
checkdata > objects are not equal; if this is expected replace base64 data with:
eNqNkl1uAzEIhK+TSHbF8GsOlPtfoTbstnlJ1RezgvH4Y7QYDx2w53g8pvGy17ShQXLqhCS/pg5XidMA
HYGOyUlUCoZod9ai88FBOFUTqwbKmi3NbCk8S4q0ktiKcnWgTW0BfWW5n4mYF5BZlPkNOhVcSpYo0B9A
4JDvgXZDMrTfRDl66L7oI9kvOu8VNVY0BOLmFrRn2LWAoUa/wIYOzHm15CaFd1C7lgKhWlzk1XdK77lF
97m4Z4TwAZw7vnpMXFY7QPtGO+0MoolTljcFqFLZy+YllStP0jKVK56ZeVbwLcl2hTZPdrDCQh2wVcD0
NfavgpDygZwKjIOCrZb3chw+FbQE+KPQe9EPp/zj5LezvZ/j+Q14zqVr

AssertionError()
  File "examples/drivencavity", line 84, in unittest
    qe0DvRUYct7B2fr328qvHwXKg10=''' )
  File "nutils/util.py", line 524, in run
    func( **kwargs )

On amd64 testing commit c207d13 (two different inequalities using python3!):

$ python3 examples/drivencavity unittest --tbexplore=False --verbose=3
checkdata > compare > #1 > #3 > non equal to 4 digits: 6.796537e-16 != 0.000000e+00
checkdata > objects are not equal; if this is expected replace base64 data with:
eNpVUuuNgzEIW6eVkhPmzUDdf4VL4Kt09yc4xBiDgvXSBXuv12sDFp9ta4NiQEnIZ/uCaieYC/0gJjIg
c6jM1k9a5DdaVcetYjzUyClGUjb1EXH1EYN7DgjXBlHgNpDStQK06teqnER3Q/pfg0C2klQ1zx4ecswo
VV5dIaUbN+QW+mIRfxL62BanZ5Dwx2bof7+GsgYez+a+Rp3uALFMqSP4Go/lqdX3rDagAbt3seGfyRkN
nEguSEInZuSrJNoEWPWONhNmGIuptbKOquajddo0EPYuRhv3lWzdNg71RsY0hd0mZ09mvR76WeeroGiy
bmd+4Hjh0xlnv32+1/sXf/SB/g==

AssertionError()
  File "examples/drivencavity", line 92, in unittest
    nf2B44XPZNjU37/5EIEz''' )
  File "nutils/util.py", line 524, in run
    func( **kwargs )

$ python3 examples/drivencavity unittest --tbexplore=False --verbose=3
checkdata > compare > #1 > non matching lenghts: 4 != 6
checkdata > compare > #1 > #3 > non equal to 4 digits: 2.182829e-15 != 0.000000e+00
checkdata > objects are not equal; if this is expected replace base64 data with:
eNpVUuuNgzEIW6eVkhPmzUDdf4VL4Kt09yc4xBiDgvXSBXuv12sDFp9ta4NiQEnIZ/uCaieYC/0gJjIg
c6jM1k9a5DdaVcetYjzUyClGUjb1EXH1EYN7DgjXBlHgNpDStQK06teqnER3Q/pfg0C2klQ1zx4ecswo
VV5dIaUbN+QW+mIRfxL62BanZ5Dwx2bof7+GsgYez+a+Rp3uALFMqSP4Go/lqdX3rDagAbt3seGfyRkN
nEguSEInZuSrJNoEWPWONhNmGIuptbKOquajddo0EPYuRhv3lWzdNg71RsY0hd0mZ09mvR76WeeroGiy
bmd+4Hjh0xln/3/O5r5/AR2Ng98=

AssertionError()
  File "examples/drivencavity", line 92, in unittest
    nf2B44XPZNjU37/5EIEz''' )
  File "nutils/util.py", line 524, in run
    func( **kwargs )

redesign the viewer

traceback explorer broken?

Somewhere between nutils/master and nutils/devel the traceback explorer seems to have broken. To be precise it happened at 170ec59

It barely seems necessary to add a script to trigger this, but here it comes

from nutils import *
def temp():
  raise ValueError
util.run( temp )

Is a unittest for this a good idea?

Remove preconditioner default

Consider to revert because diag does not always work, e.g. for mixed systems (stokes etc)

Possible bug Makefile nutils

When I try to compile the nutils _numeric,c file I encountered the following error in inside the folder nutils/nutils

$ make
python setup.py build
Traceback (most recent call last):
  File "/usr/lib/python2.7/site.py", line 70, in <module>
    import traceback
  File "traceback.py", line 1, in <module>
    from . import core
ValueError: Attempted relative import in non-package
make: *** [_numeric.so] Error 1

It runs when I replace from . import core with import core
I don't know if this is a proper fix for the problem but after that it worked again.

Best,
Rody

Hard disk space is used when using shared zeros (nprocs>=1)

delayed integration object

When integrating an integrand with placeholder arguments return a DelayedIntegral object, name up for discussion, instead of evaluating the integral immediately. DelayedIntegral objects can be summed if the shapes match. A DelayedIntegral object can be evaluated via e.g. .eval, which requires values for all placeholder arguments.

test issue

Just trying it out.

Error reported

I got this error while executing make test

python setup.py build clean
running build
running build_py
copying nutils/mesh.py -> build/lib.macosx-10.5-x86_64-2.7/nutils
running build_ext
running clean
cp -vu build/lib_/nutils/numeric.so nutils
cp: illegal option -- u
usage: cp [-R [-H | -L | -P]] [-fi | -n] [-apvX] source_file target_file
cp [-R [-H | -L | -P]] [-fi | -n] [-apvX] source_file ... target_directory
make: ** [dev] Error 64

Bug in integration over trimmed boundary

We get the following error:

nutils v1.dev.3ab03374

curved main \
  --R=2.0 \
  --r=1.0 \
  --l=1.0 \
  --lmbda=1.0 \
  --mu=1.0 \
  --local=False \
  --mmr=0 \
  --degree=1 \
  --solvetol=1e-10 \
  --h_v=[ 2.   1.   0.5]

start Wed Sep 16 13:34:59 2015

starting with h_n = 0
trim · elem 43/100 (42%)
trim · elem 65/100 (64%)
created geometry in 4.34107112885
created volumetric operators in 0.275549888611
finished volume integration in 0.912701129913
created boundary operators in 1.39767408371
ValueError('tuple.index(x): x not in tuple',)
  File "../nutils/nutils/topology.py", line 1311, in boundary
    iedge = ref.edge_transforms.index(tail)
  File "../nutils/nutils/cache.py", line 26, in property_getter
    value = f( self )
  File "../nutils/nutils/topology.py", line 1303, in boundary
    basebtopo = self.basetopo.boundary
  File "../nutils/nutils/cache.py", line 26, in property_getter
    value = f( self )
  File "curved", line 309, in main
    stability_tang_matrix_bottom = domain.boundary['trim_bottom'].integrate( [ flux_norm_matrix,
  File "../nutils/nutils/util.py", line 525, in run
    func( **kwargs )

When running the code below:

#!/usr/bin/env python

from nutils import *
import scipy.sparse.linalg
import time

def precon(M,tol_orth=.99,tol_rem=1e-12):
  islands = []
  p = numpy.ones(M.shape[0],dtype=bool)
  D = scipy.sparse.lil_matrix(scipy.sparse.diags(1./numpy.sqrt(M.diagonal()),0))
  i_v,j_v,dummy = scipy.sparse.find( numpy.abs( scipy.sparse.tril( D.dot(M.dot(D.T)) - scipy.sparse.diags( D.dot(M.dot(D.T)).diagonal(),0 ) ) ) > tol_orth )
  while len(i_v): # while the matrix requires orthonormalisation
    ## sort indices
    for index in zip( i_v, j_v ):
      index = set( index )
      for island in list(islands):
        if island & index:
          islands.remove( island )
          index |= island
      islands.append( index )
    ## orthonormalise
    for island in islands:
      island_unsort = list( island )
      other_indices = list( set(numpy.arange(M.shape[0])) - island )
      intersections = numpy.ravel((numpy.abs(D.dot(M.dot(D.T))[:,other_indices][island_unsort,:])>tol_rem).sum(1))
      island = []
      for island_index in intersections.argsort():
        island.append( island_unsort[island_index] )
      M_p = M[:,island][island,:]
      D_p = numpy.diagflat(1./numpy.sqrt(M_p.diagonal()),0)
      for ind_1 in numpy.arange( len( island ) ):
        for ind_2 in numpy.arange( ind_1 ):
          D_p[ind_1,:] -= D_p[ind_2,:]*D_p[ind_2,:].dot(M_p.dot(D_p[ind_1,:]))
        if D_p[ind_1,:].dot(M_p.dot(D_p[ind_1,:])) <= tol_rem:
          D_p[ind_1,:] = 0
          p[island[ind_1]] = False
        else:
          D_p[ind_1,:] = D_p[ind_1,:]/numpy.sqrt(D_p[ind_1,:].dot(M_p.dot(D_p[ind_1,:])))
      D[:,island][island,:] = D_p
    i_v,j_v,dummy = scipy.sparse.find( numpy.abs( scipy.sparse.tril( D.dot(M.dot(D.T)) - scipy.sparse.diags( D.dot(M.dot(D.T)).diagonal(),0 ) ) ) > tol_orth )
  return D[p,:]

def plotmesh ( domain, geom, disp=0., norm_stress=0., mag_stress=0., name='curved' ):
  points, disp, norm_stress, mag_stress = domain.boundary.simplex.elem_eval( [geom,disp,norm_stress,mag_stress], ischeme='vtk', separate=True )
  with plot.VTKFile(name+'_boundary') as vtkfile:
    vtkfile.unstructuredgrid( points )
    vtkfile.pointdataarray( 'disp', disp )
    vtkfile.pointdataarray( 'norm_stress', norm_stress )
    vtkfile.pointdataarray( 'mag_stress', mag_stress )
#  points, displacement, variable = domain         .simplex.elem_eval( [geom,disp,var], ischeme='vtk', separate=True )
#  with plot.VTKFile(name+'_volume') as vtkfile:
#    vtkfile.unstructuredgrid( points )
#    vtkfile.pointdataarray( 'disp', disp )
#    vtkfile.pointdataarray( 'norm_stress', norm_stress )
#    vtkfile.pointdataarray( 'mag_stress', mag_stress )
  return

def main( R     =  2. , #Radius curve
          r     =  1. , #Radius rod
          l     =  1. , #Length rod
          lmbda =  1. , #First Lame parameter
          mu    =  1. , #Second Lame parameter
          local    = False  , #Stabilisation
          mmr      = 0      , #Minmaxrefine
          degree   = 1      , #Degree
          solvetol = 1e-10  , #Solvetol
          h_v      = numpy.array([2.,1.,.5]) #Gridsize
        ):

  stress = library.Hooke(lmbda=lmbda,mu=mu)

  energy_norm_v = numpy.zeros(h_v.shape)
  energy_tang_v = numpy.zeros(h_v.shape)

  for h_n in numpy.arange(h_v.shape[0]):
    h = h_v[h_n]
    log.user('starting with h_n =',h_n)
    t_old = time.time()

    ni = numpy.ceil(  l  /h)  #Number of inner elements 
    no = numpy.ceil((R+r)/h)  #Number of outer elements 
    nz = numpy.ceil(  r  /h)  #Number of z elements

    domain, geom = mesh.rectilinear([ numpy.linspace( -h*(ni+degree), h*(no+degree), ni+no+2*degree+1 ),
                                      numpy.linspace( -h*(ni+degree), h*(no+degree), ni+no+2*degree+1 ),
                                      numpy.linspace( -h*(nz+degree), h*(nz+degree),  2*nz+2*degree+1 ) ])

    x, y, z = geom

    bc_norm = function.asarray([ 1.,  0., 0.   ])
    bc_tang = function.asarray([ 0., -z , y-R  ])

    dxy = function.piecewise( function.arctan2(y,x), [0,numpy.pi/2], (x-R)**2                                                ,
                                                                     x**2 + y**2 + R**2 - 2*R*function.sqrt(x**2+y**2),
                                                                     (y-R)**2                                                  )
    lvl_rad    = r**2 - dxy - z**2
    lvl_bottom = l + y
    lvl_left   = l + x
    domain = domain.trim( lvl_rad   , maxrefine = mmr+h_v.shape[0]-h_n-1, name = 'trim_rad'    )
    domain = domain.trim( lvl_bottom, maxrefine = mmr+h_v.shape[0]-h_n-1, name = 'trim_bottom' )
    domain = domain.trim( lvl_left  , maxrefine = mmr+h_v.shape[0]-h_n-1, name = 'trim_left'   )

    tang_tensor    = function.eye( 3 ) - geom.normal()[:,_]*geom.normal()[_,:]
    funcsp         = domain.basis( 'spline', degree=degree ).vector( 3 )
    symgrad_funcsp = funcsp.symgrad(geom)
    stresssp       = stress(symgrad_funcsp)
    log.user('created geometry in',time.time()-t_old)

    ##############################
    # Volume integration         #
    ##############################

    t_old = time.time()
    elasticity_matrix = function.outer( symgrad_funcsp, stresssp ).sum([-1,-2])
    log.user('created volumetric operators in',time.time()-t_old)

    t_old = time.time()
    volume_matrix, locs, vols = domain.integrate( [ elasticity_matrix,
                                                    (symgrad_funcsp**2).sum([-1,-2])[:,_]*geom[_,:],
                                                    (symgrad_funcsp**2).sum([-1,-2])
                                                  ], geometry=geom, ischeme='gauss'+str(2*degree) )
    volume_matrix = volume_matrix.toscipy()
    locs          = locs.toarray()
    xlocs = locs[2*funcsp.shape[0]/3::                  ,0]/vols[2*funcsp.shape[0]/3::                  ]
    ylocs = locs[                   :  funcsp.shape[0]/3,1]/vols[                   :  funcsp.shape[0]/3]
    zlocs = locs[  funcsp.shape[0]/3:2*funcsp.shape[0]/3,2]/vols[  funcsp.shape[0]/3:2*funcsp.shape[0]/3]
    xfacts = (ylocs-numpy.average(ylocs)) * volume_matrix.diagonal()[                   :  funcsp.shape[0]/3]
    yfacts = (zlocs-numpy.average(zlocs)) * volume_matrix.diagonal()[  funcsp.shape[0]/3:2*funcsp.shape[0]/3]
    zfacts = (xlocs-numpy.average(xlocs)) * volume_matrix.diagonal()[2*funcsp.shape[0]/3::                  ]
    x_index  =                       numpy.argmax(volume_matrix.diagonal()[                    :   funcsp.shape[0]/3])
    y_index  =   funcsp.shape[0]/3 + numpy.argmax(volume_matrix.diagonal()[  funcsp.shape[0]/3 : 2*funcsp.shape[0]/3])
    z_index  = 2*funcsp.shape[0]/3 + numpy.argmax(volume_matrix.diagonal()[2*funcsp.shape[0]/3 : :                  ])
    x_index1 =                       numpy.argmin(xfacts); x_index2 =                       numpy.argmax(xfacts)
    y_index1 =   funcsp.shape[0]/3 + numpy.argmin(yfacts); y_index2 =   funcsp.shape[0]/3 + numpy.argmax(yfacts)
    z_index1 = 2*funcsp.shape[0]/3 + numpy.argmin(zfacts); z_index2 = 2*funcsp.shape[0]/3 + numpy.argmax(zfacts)
    log.user('finished volume integration in',time.time()-t_old)

    ##############################
    # Setting stability constants#
    ##############################

    if local:

      if h_n == 0: 

        #Create local bases indices
        x_power = numpy.reshape(numpy.tile(numpy.arange(degree+1),(degree+1,degree+1,1)).transpose((0,1,2)),-1)
        y_power = numpy.reshape(numpy.tile(numpy.arange(degree+1),(degree+1,degree+1,1)).transpose((1,2,0)),-1)
        z_power = numpy.reshape(numpy.tile(numpy.arange(degree+1),(degree+1,degree+1,1)).transpose((2,0,1)),-1)
        basis_indices_mu = []; basis_indices_lmbda = []
        for i in numpy.arange((degree+1)**3):
          basis_indices_mu   .append((x_power[i],y_power[i],z_power[i]))
          basis_indices_lmbda.append((x_power[i],y_power[i],z_power[i]))
        basis_indices_mu   .remove(( 0     , 0     , 0      ))
        basis_indices_mu   .remove(( 1     , 0     , 0      ))
        basis_indices_mu   .remove(( 0     , 1     , 0      ))
        basis_indices_mu   .remove(( 0     , 0     , 1      ))
        basis_indices_lmbda.remove(( degree, degree, degree ))

        #Initialise
        beta_mu_norm_dict = {}; beta_mu_tang_dict = {}; beta_lmbda_dict = {}

        #Index Dirichlet boundary elements
        elem_belems = {}
        for belem in domain.boundary['trim_bottom,trim_left']:
          etrans = belem.transform.lookup( domain.edict )
          ielem  = domain.edict[ etrans ]
          elem_belems.setdefault(ielem,[]).append( belem )

        #Loop over Dirichlet boundary elements
        for ielem, belems in log.iter( 'boundary element', elem_belems.items() ):
          velem = domain.elements[ielem]

          #Create local topology
          vdomain = topology.Topology( [velem] )
          bdomain = topology.Topology( belems )

          #Calculate centre of mass
          volume, geom_sep = vdomain.integrate( [ 1., geom ], geometry=geom, ischeme='gauss1' )
          geom_local = geom - geom_sep / volume

          #Create local bases
          basis_blocks = [ f * function.eye(3) for f in function.product( function.power( geom_local[_,:], basis_indices_mu ), axis=-1 ) ]
          basis_blocks.append( function.diagonalize(geom_local) )
          basis_blocks.append([ [ geom_local[1], geom_local[0], 0             ],
                                [ 0            , geom_local[2], geom_local[1] ],
                                [ geom_local[2], 0            , geom_local[0] ] ])
          basis_stab_mu = function.concatenate( basis_blocks, axis=0 )

          basis_stab_lmbda = function.product( function.power( geom_local[_,:], basis_indices_lmbda ), axis=-1 )

          #Operators
          symgrad_basis_mu   = basis_stab_mu.symgrad(geom)
          E_mu_norm_operator = function.outer(  symgrad_basis_mu.dotnorm(geom).dotnorm(geom)                        )
          E_mu_tang_operator = function.outer( (symgrad_basis_mu.dotnorm(geom)[:,_,:] * tang_tensor[_,:,:]).sum(-1) ).sum(-1)
          E_lmbda_operator   = function.outer( basis_stab_lmbda )
          A_mu_operator      = function.outer( symgrad_basis_mu ).sum([-1,-2])
          A_lmbda_operator   = E_lmbda_operator

          #Integrate
          E_mu_norm_elem, E_mu_tang_elem, E_lmbda_elem = bdomain.integrate( [ E_mu_norm_operator, E_mu_tang_operator, E_lmbda_operator ], geometry=geom, ischeme='gauss'+str(2*degree) )
          A_mu_elem, A_lmbda_elem                      = vdomain.integrate( [ A_mu_operator, A_lmbda_operator                          ], geometry=geom, ischeme='gauss'+str(2*degree) )
          E_mu_norm_elem = E_mu_norm_elem.toscipy(); E_mu_tang_elem = E_mu_tang_elem.toscipy(); E_lmbda_elem = E_lmbda_elem.toscipy()
          A_mu_elem = A_mu_elem.toscipy(); A_lmbda_elem = A_lmbda_elem.toscipy()

          #Eigenvalue problem
#          D = precon(A_mu_elem   ); D = D.toarray()
#          beta_mu_norm_dict[velem.transform] = 4*mu   *max(scipy.linalg.eigh(D.dot(E_mu_norm_elem.dot(D.T)),D.dot(A_mu_elem   .dot(D.T)),eigvals_only=True))
#          beta_mu_tang_dict[velem.transform] = 4*mu   *max(scipy.linalg.eigh(D.dot(E_mu_tang_elem.dot(D.T)),D.dot(A_mu_elem   .dot(D.T)),eigvals_only=True))
#          D = precon(A_lmbda_elem); D = D.toarray()
#          beta_lmbda_dict  [velem.transform] = 2*lmbda*max(scipy.linalg.eigh(D.dot(E_lmbda_elem  .dot(D.T)),D.dot(A_lmbda_elem.dot(D.T)),eigvals_only=True))
          D = precon(A_mu_elem   ); D = D.toarray()
          beta_h_mu_norm = 4*mu   *max(scipy.linalg.eigh(D.dot(E_mu_norm_elem.dot(D.T)),D.dot(A_mu_elem   .dot(D.T)),eigvals_only=True))*h
          beta_h_mu_tang = 4*mu   *max(scipy.linalg.eigh(D.dot(E_mu_tang_elem.dot(D.T)),D.dot(A_mu_elem   .dot(D.T)),eigvals_only=True))*h
          D = precon(A_lmbda_elem); D = D.toarray()
          beta_h_lmbda   = 2*lmbda*max(scipy.linalg.eigh(D.dot(E_lmbda_elem  .dot(D.T)),D.dot(A_lmbda_elem.dot(D.T)),eigvals_only=True))*h

#          D = scipy.sparse.csr_matrix(precon(A_mu_elem   ))
#          (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_mu_norm_elem.dot(D.T)), k=1, M=D.dot(A_mu_elem   .dot(D.T)), which='LM' )
#          beta_mu_norm_dict[velem.transform] = 4*mu    * eigenvalue
#          (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_mu_tang_elem.dot(D.T)), k=1, M=D.dot(A_mu_elem   .dot(D.T)), which='LM' )
#          beta_mu_tang_dict[velem.transform] = 4*mu    * eigenvalue
#          D = scipy.sparse.csr_matrix(precon(A_lmbda_elem))
#          (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_lmbda_elem  .dot(D.T)), k=1, M=D.dot(A_lmbda_elem.dot(D.T)), which='LM' )
#          beta_lmbda_dict  [velem.transform] = 2*lmbda * eigenvalue
#          D = scipy.sparse.csr_matrix(precon(A_mu_elem   ))
#          (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_mu_norm_elem.dot(D.T)), k=1, M=D.dot(A_mu_elem   .dot(D.T)), which='LM' )
#          beta_h_mu_norm = 4*mu    * eigenvalue * h
#          (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_mu_tang_elem.dot(D.T)), k=1, M=D.dot(A_mu_elem   .dot(D.T)), which='LM' )
#          beta_h_mu_tang = 4*mu    * eigenvalue * h
#          D = scipy.sparse.csr_matrix(precon(A_lmbda_elem))
#          (eigenvalue,),dummy = scipy.sparse.linalg.eigsh( A=D.dot(E_lmbda_elem  .dot(D.T)), k=1, M=D.dot(A_lmbda_elem.dot(D.T)), which='LM' )
#          beta_h_lmbda   = 2*lmbda * eigenvalue * h

          break 

#        #Assign 
#        beta_mu_norm = function.Elemwise( beta_mu_norm_dict, () )
#        beta_mu_tang = function.Elemwise( beta_mu_tang_dict, () )
#        beta_lmbda   = function.Elemwise( beta_lmbda_dict  , () )

      beta_mu_norm = beta_h_mu_norm/h
      beta_mu_tang = beta_h_mu_tang/h
      beta_lmbda   = beta_h_lmbda  /h

      ##############################
      # Boundary integration       #
      ##############################

      t_old = time.time()
      #Operators
      flux_norm_matrix    = -function.outer( stresssp.dotnorm(geom).dotnorm(geom) , funcsp.dotnorm(geom) )
      flux_norm_matrix    = flux_norm_matrix + flux_norm_matrix.T
      flux_norm_vector    = -stresssp.dotnorm(geom).dotnorm(geom) * bc_norm.dotnorm(geom)
      penalty_norm_matrix = (beta_mu_norm+beta_lmbda) * function.outer( funcsp.dotnorm(geom) )
      penalty_norm_vector = (beta_mu_norm+beta_lmbda) * funcsp.dotnorm(geom) * bc_norm.dotnorm(geom)
      flux_tang_matrix    = -function.outer( tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,:,_], funcsp[:,_,:] ).sum([-1,-2])
      flux_tang_matrix    = flux_tang_matrix + flux_tang_matrix.T
      flux_tang_vector    = -(tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,:,_] * bc_tang[_,_,:] ).sum([-1,-2])
      penalty_tang_matrix = beta_mu_tang * function.outer( ( tang_tensor[_,:,:] * funcsp[:,_,:] ).sum(-1) ).sum(-1)
      penalty_tang_vector = beta_mu_tang * ( tang_tensor[_,:,:] * funcsp[:,:,_] * bc_tang[_,_,:] ).sum([-1,-2])
      log.user('created boundary operators in',time.time()-t_old)

      t_old = time.time()
      #Integrate
      boundary_norm_matrix_bottom,                     boundary_tang_matrix_bottom                     = domain.boundary['trim_bottom'].integrate( [ flux_norm_matrix + penalty_norm_matrix,
                                                                                                                                                     flux_tang_matrix + penalty_tang_matrix,
                                                                                                                                                   ], geometry=geom, ischeme='gauss'+str(2*degree) )

      log.user('finished bottom integration in',time.time()-t_old)

      t_old = time.time()
      boundary_norm_matrix_left, boundary_norm_vector, boundary_tang_matrix_left, boundary_tang_vector = domain.boundary['trim_left'  ].integrate( [ flux_norm_matrix + penalty_norm_matrix,
                                                                                                                                                     flux_norm_vector + penalty_norm_vector,
                                                                                                                                                     flux_tang_matrix + penalty_tang_matrix,
                                                                                                                                                     flux_tang_vector + penalty_tang_vector
                                                                                                                                                   ], geometry=geom, ischeme='gauss'+str(2*degree) )

      log.user('finished left integration in',time.time()-t_old)

      boundary_norm_matrix = boundary_norm_matrix_bottom.toscipy() + boundary_norm_matrix_left.toscipy() + boundary_tang_matrix_bottom.toscipy()
      boundary_tang_matrix = boundary_tang_matrix_bottom.toscipy() + boundary_tang_matrix_left.toscipy()

    else:

      t_old = time.time()
      #Operators 
      flux_norm_matrix       = -function.outer( stresssp.dotnorm(geom).dotnorm(geom) , funcsp.dotnorm(geom) )
      flux_norm_matrix       = flux_norm_matrix + flux_norm_matrix.T
      flux_norm_vector       = -stresssp.dotnorm(geom).dotnorm(geom) * bc_norm.dotnorm(geom)
      penalty_norm_matrix    = function.outer( funcsp.dotnorm(geom) )
      penalty_norm_vector    = funcsp.dotnorm(geom) * bc_norm.dotnorm(geom)
      stability_norm_matrix  = function.outer( stresssp.dotnorm(geom).dotnorm(geom) )
      flux_tang_matrix       = -function.outer( tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,:,_], funcsp[:,_,:] ).sum([-1,-2])
      flux_tang_matrix       = flux_tang_matrix + flux_tang_matrix.T
      flux_tang_vector       = -(tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,:,_] * bc_tang[_,_,:] ).sum([-1,-2])
      penalty_tang_matrix    = function.outer( ( tang_tensor[_,:,:] * funcsp[:,_,:] ).sum(-1) ).sum(-1)
      penalty_tang_vector    = ( tang_tensor[_,:,:] * funcsp[:,:,_] * bc_tang[_,_,:] ).sum([-1,-2])
      stability_tang_matrix  = function.outer( ( tang_tensor[_,:,:] * stresssp.dotnorm(geom)[:,_,:] ).sum(-1) ).sum(-1)
      log.user('created boundary operators in',time.time()-t_old)

      t_old = time.time()
      #Integrate
      flux_norm_matrix_bottom, \
      penalty_norm_matrix_bottom, \
      stability_norm_matrix_bottom, \
      flux_tang_matrix_bottom, \
      penalty_tang_matrix_bottom, \
      stability_tang_matrix_bottom = domain.boundary['trim_bottom'].integrate( [ flux_norm_matrix, 
                                                                                 penalty_norm_matrix,
                                                                                 stability_norm_matrix,
                                                                                 flux_tang_matrix,
                                                                                 penalty_tang_matrix,
                                                                                 stability_tang_matrix    
                                                                               ], geometry=geom, ischeme='gauss'+str(2*degree) )
      log.user('finished bottom integration in',time.time()-t_old)

      t_old = time.time()
      flux_norm_matrix_left, \
      flux_norm_vector, \
      penalty_norm_matrix_left, \
      penalty_norm_vector, \
      stability_norm_matrix_left, \
      flux_tang_matrix_left, \
      flux_tang_vector, \
      penalty_tang_matrix_left, \
      penalty_tang_vector, \
      stability_tang_matrix_left   = domain.boundary['trim_left'  ].integrate( [ flux_norm_matrix,
                                                                                 flux_norm_vector,
                                                                                 penalty_norm_matrix,
                                                                                 penalty_norm_vector,
                                                                                 stability_norm_matrix,
                                                                                 flux_tang_matrix,
                                                                                 flux_tang_vector,
                                                                                 penalty_tang_matrix,
                                                                                 penalty_tang_vector,
                                                                                 stability_tang_matrix    
                                                                               ], geometry=geom, ischeme='gauss'+str(2*degree) )
      log.user('finished left integration in',time.time()-t_old)

      flux_norm_matrix      = flux_norm_matrix_bottom     .toscipy() + flux_norm_matrix_left     .toscipy() + flux_tang_matrix_bottom     .toscipy()
      penalty_norm_matrix   = penalty_norm_matrix_bottom  .toscipy() + penalty_norm_matrix_left  .toscipy() + penalty_tang_matrix_bottom  .toscipy()
      stability_norm_matrix = stability_norm_matrix_bottom.toscipy() + stability_norm_matrix_left.toscipy() + stability_tang_matrix_bottom.toscipy() 
      flux_tang_matrix      = flux_tang_matrix_bottom     .toscipy() + flux_tang_matrix_left     .toscipy()
      penalty_tang_matrix   = penalty_tang_matrix_bottom  .toscipy() + penalty_tang_matrix_left  .toscipy()
      stability_tang_matrix = stability_tang_matrix_bottom.toscipy() + stability_tang_matrix_left.toscipy()

      #Remove singularities
      p = numpy.ones(funcsp.shape[0],dtype=bool)
      p[ x_index1 ] = False
      p[ x_index2 ] = False
      p[ y_index1 ] = False
      p[ y_index2 ] = False
      p[ z_index1 ] = False
      p[ z_index2 ] = False

      t_old = time.time()
      #Eigenvalue problem
      A_matrix      = volume_matrix        [:,p][p,:]
      E_norm_matrix = stability_norm_matrix[:,p][p,:]
      E_tang_matrix = stability_tang_matrix[:,p][p,:]
      D = precon(A_matrix).toarray()
      log.user('assembled eigenvalue problems in',time.time()-t_old)

      t_old = time.time()
      beta_norm = 2*max(scipy.linalg.eigh(D.dot(E_norm_matrix.dot(D.T)),D.dot(A_matrix.dot(D.T)),eigvals_only=True))
      beta_tang = 2*max(scipy.linalg.eigh(D.dot(E_tang_matrix.dot(D.T)),D.dot(A_matrix.dot(D.T)),eigvals_only=True))
#      (beta_norm,),dummy = 2*scipy.sparse.linalg.eigsh( A=D.dot(E_norm_matrix.dot(D.T)), k=1, M=D.dot(A_matrix.dot(D.T)), which='LM' )
#      (beta_tang,),dummy = 2*scipy.sparse.linalg.eigsh( A=D.dot(E_tang_matrix.dot(D.T)), k=1, M=D.dot(A_matrix.dot(D.T)), which='LM' )
      log.user('solved eigenvalue problems in',time.time()-t_old)

      #Assemble
      boundary_norm_matrix = flux_norm_matrix + beta_norm * penalty_norm_matrix
      boundary_norm_vector = flux_norm_vector + beta_norm * penalty_norm_vector
      boundary_tang_matrix = flux_tang_matrix + beta_tang * penalty_tang_matrix
      boundary_tang_vector = flux_tang_vector + beta_tang * penalty_tang_vector

    ##############################
    # Solve                      #
    ##############################

    #Remove singularities
    p_norm   = numpy.ones(funcsp.shape[0],dtype=bool)
    p_tang   = numpy.ones(funcsp.shape[0],dtype=bool)
    p_tang[ x_index2 ] = False

    t_old = time.time() 
    #Assemble
    full_norm_matrix = (volume_matrix + boundary_norm_matrix)[:,p_norm][p_norm,:]
    full_norm_vector = boundary_norm_vector[p_norm]
    full_tang_matrix = (volume_matrix + boundary_tang_matrix)[:,p_tang][p_tang,:]
    full_tang_vector = boundary_tang_vector[p_tang]
    log.user('assembled full systems in',time.time()-t_old)

    #Solve
    t_old = time.time()
    D = scipy.sparse.csr_matrix(precon(full_norm_matrix))
    lhs_free = D.T.dot(scipy.sparse.linalg.cg(D.dot(full_norm_matrix.dot(D.T)),D.dot(full_norm_vector),tol=solvetol,maxiter=10000)[0])
    lhs_norm = numpy.zeros(funcsp.shape[0]); lhs_norm[p_norm] = lhs_free
    log.user('solved norm in',time.time()-t_old)
    t_old = time.time()
    D = scipy.sparse.csr_matrix(precon(full_tang_matrix))
    lhs_free = D.T.dot(scipy.sparse.linalg.cg(D.dot(full_tang_matrix.dot(D.T)),D.dot(full_tang_vector),tol=solvetol,maxiter=10000)[0])
    lhs_tang = numpy.zeros(funcsp.shape[0]); lhs_tang[p_tang] = lhs_free
    log.user('solved tang in',time.time()-t_old)

    ##############################
    # Postprocessing             #
    ##############################

    #Convergence check
    energy_norm_v[h_n] = lhs_norm.dot(volume_matrix.dot(lhs_norm))
    energy_tang_v[h_n] = lhs_tang.dot(volume_matrix.dot(lhs_tang))
    if h_n>0:
      with plot.PyPlot('energy_normal_' ) as plt:
        plt.loglog(h_v[:h_n],numpy.abs(energy_norm_v[:h_n]-energy_norm_v[1:h_n+1]),c='r',marker='o',markersize=10,markeredgecolor='r',linewidth=0,fillstyle='full')
        plt.loglog(h_v[:h_n],numpy.abs(energy_norm_v[h_n-1]-energy_norm_v[h_n])/h_v[h_n-1]**(2*degree)*h_v[:h_n]**(2*degree),c='k',linewidth=1,)
        plt.grid(True)
        #plt.legend(fontsize=20)
        plt.xlabel(r'$h$',fontsize=20)
        plt.ylabel(r'$E(v^h)$',fontsize=20)
      log.user('normal' ,numpy.abs(energy_norm_v[:h_n]-energy_norm_v[1:h_n+1]))
      with plot.PyPlot('energy_tangent_') as plt:
        plt.loglog(h_v[:h_n],numpy.abs(energy_tang_v[:h_n]-energy_tang_v[1:h_n+1]),c='r',marker='o',markersize=10,markeredgecolor='r',linewidth=0,fillstyle='full')
        plt.loglog(h_v[:h_n],numpy.abs(energy_tang_v[h_n-1]-energy_tang_v[h_n])/h_v[h_n-1]**(2*degree)*h_v[:h_n]**(2*degree),c='k',linewidth=1,)
        plt.grid(True)
        #plt.legend(fontsize=20)
        plt.xlabel(r'$h$',fontsize=20)
        plt.ylabel(r'$E(v^h)$',fontsize=20)
      log.user('tangent',numpy.abs(energy_tang_v[:h_n]-energy_tang_v[1:h_n+1]))

    #Plotting
    plotmesh( domain, geom, funcsp.dot(lhs_norm), stresssp.dot(lhs_norm).dotnorm(geom), (stresssp.dot(lhs_norm)**2).sum([-1,-2]), 'h_'+str(int(h_n))+'_normal'  )
    plotmesh( domain, geom, funcsp.dot(lhs_tang), stresssp.dot(lhs_tang).dotnorm(geom), (stresssp.dot(lhs_tang)**2).sum([-1,-2]), 'h_'+str(int(h_n))+'_tangent' )

  return

util.run( main )

Changed numpy to function

DG function space with degree = 1

I have an issue when trying to run a simulation when I write domian.discontfunc(degree=1) on Finity (I did not make the switch to Nutils yet). I get the following error:

IndexError: index 1 is out of bounds for axis 1 with size 1
File "../Finity/finity/element.py", line 1609, in bernstein_poly
root = ( poly[k,k+1] * (k_2-n+1) ) / (k+1)
File "../Finity/finity/topology.py", line 837, in
stdelem = util.product( element.PolyLine( element.PolyLine.bernstein_poly( d ) ) for d in degree )
File "../Finity/finity/util.py", line 186, in product
prod = seq.next()
File "../Finity/finity/topology.py", line 837, in discontfunc
stdelem = util.product( element.PolyLine( element.PolyLine.bernstein_poly( d ) ) for d in degree )
File "../Finity/finity/core.py", line 24, in wrapped
value = func( *args )
File "Euler1DExplicit", line 40, in main
Rspace = domain.discontfunc(degree = degree)
File "../Finity/finity/util.py", line 474, in run
func( *_kwargs )

reusing arrayfunc structure

Nutils is currently unable to reuse function trees such as occurring in Newton/Picard iterations or time stepping, where integrands have identical structure but coefficients are changing. Requires some thinking but it's an important potential for optimization.

Evaluation of trimmed boundaries

When a domain is trimmed, is seems that only the boundary created in the last trimming operation is evaluated as the trimmed domain's boundary. Both plotting the boundary as evaluating the length of the boundary yield erroneous results. I have tried this with a unit square that is rotated around its center. An exception to these errors is the unrotated square, which has all the boundaries, but also not with the correct length of the separate parts.

#! /usr/bin/env python3

from nutils import *

def main( nelems = 8,
          nsteps = 101,
        ):

  theta_v = numpy.pi/4*numpy.linspace(0,1,nsteps) #angles of rotation

  verts = numpy.linspace( -1, 1, 2*nelems+1 )
  untrimmed_domain, geom = mesh.rectilinear([ verts,verts ])

  for theta_n in numpy.arange(theta_v.shape[0]):
    log.user('__________________________________________________')
    log.user('##          computing theta',theta_n,'of',nsteps)
    log.user('__________________________________________________')
    theta = theta_v[theta_n]

    Q = function.asarray([[ numpy.cos(theta),-numpy.sin(theta) ],
                          [ numpy.sin(theta), numpy.cos(theta) ]]) #rotation tensor
    geom_rel = (Q.T*geom[_,:]).sum(-1) #coordinates relative to rotated object

    #trim domain w.r.t. the boundaries of the rotated object
    domain = untrimmed_domain.trim( .5+geom_rel[0], maxrefine=3, ndivisions=15, name='trimleft'   )
    domain =           domain.trim( .5-geom_rel[0], maxrefine=3, ndivisions=15, name='trimright'  )
    domain =           domain.trim( .5+geom_rel[1], maxrefine=3, ndivisions=15, name='trimbottom' )
    domain =           domain.trim( .5-geom_rel[1], maxrefine=3, ndivisions=15, name='trimtop'    )

    boundary_length = domain.boundary['trimleft'].integrate( 1, geometry=geom, ischeme='gauss2' ) #length of the boundary, should equal 1
    log.user('boundary length left =',boundary_length)
    boundary_length = domain.boundary['trimright'].integrate( 1, geometry=geom, ischeme='gauss2' ) #length of the boundary, should equal 1
    log.user('boundary length right =',boundary_length)
    boundary_length = domain.boundary['trimbottom'].integrate( 1, geometry=geom, ischeme='gauss2' ) #length of the boundary, should equal 1
    log.user('boundary length bottom =',boundary_length)
    boundary_length = domain.boundary['trimtop'].integrate( 1, geometry=geom, ischeme='gauss2' ) #length of the boundary, should equal 1
    log.user('boundary length top =',boundary_length)
    boundary_length = domain.boundary.integrate( 1, geometry=geom, ischeme='gauss2' ) #length of the boundary, should equal 4 for unit square
    log.user('boundary length =',boundary_length)

    #plot boundary
    x_coords, y_coords, x_rel_coords, y_rel_coords = domain.boundary.elem_eval( [geom[0], geom[1], geom_rel[0], geom_rel[1]], ischeme='bezier3', separate=False )
    with plot.PyPlot( 'abs_boundary'+str(int(theta_n)), ndigits=0 ) as plt:
      plt.plot( x_coords, y_coords )
    with plot.PyPlot( 'rel_boundary'+str(int(theta_n)), ndigits=0 ) as plt:
      plt.plot( x_rel_coords, y_rel_coords )

    #plot volume
    coords, rel_coords, vals = domain.simplex.elem_eval( [geom, geom_rel, 1], ischeme='bezier3', separate=True )
    with plot.PyPlot( 'abs_volume'+str(int(theta_n)), ndigits=0 ) as plt:
      plt.mesh( coords, vals, triangulate = 'bezier' )
    with plot.PyPlot( 'rel_volume'+str(int(theta_n)), ndigits=0 ) as plt:
      plt.mesh( rel_coords, vals, triangulate = 'bezier' )

  return

util.run(main)

add JIT backend (alongside numpy)

util.withrepr

Seems like util.withrepr gives an error when wrapped around a function without default arguments.

Html log

Run this script with verbose=4 and behold:

1. The value of the obj parameter doesn't show in the log
2. Running with user level verbosity still gives debug level output in the log

road to nutils 1.0

First: please read this blog post about upcoming changes in our branching model.

Preparing for the release of nutils 1.0, today we call on all users to check their applications against the current master branch an report on any problems they find. Please use and monitor github's issue tracker to avoid double findings. Please keep an eye on the present issue for all communication leading up to the eventual release, expected some weeks from now.

Gmesh-topology interfaces broken

The gmesh reader directly inserts an interfaces topology, bypassing the code in UnstructuredTopology.interfaces. Because of this the local coordinate systems of the both sides might not be in agreement. This means that DG type calculations on these topologies are currently broken.

Curl operator

Would it be possible to have the curl operator added to the function module?

add multipatch topologies

Documentation/working document

There should be a document that can be used to get people started with nutils. Personally I do not care about the exact form. A (elaborate) presentation would also work presumably. Of course, at some point in time, source code documentation should also be improved, but in my opinion that is not the top priority at this moment.

I'm curious to hear what others think about this. It would be great if we could come up with a concrete plan to develop this working document.

Of course the BDFL will play an essential role in the development of this document, but it is not solely his responsibility. All nutils users should be committed to work on this.

allow variable number of external arguments in serialized function tree

As serialization (and taking derivatives and optimizing the function tree) can be time consuming we want to reuse a serialized function tree whenever possible. To enable this, we should generalize the list of hardcoded external arguments or tokens of a serialized function tree to a variable number of arguments. When serializing a function all external arguments will be replaced by placeholders. Calling the serialized function (via eval) requires specifying a value for each placeholder.

Bug in integration of single trimmed element boundary

When a single trimmed element's boundary is integrated, the outcome is incorrect for small elements. Plotting the boundary and volume of the trimmed element seems to uncover what's happening. I think the figures generated by the attached file will make the problem very clear, otherwise call me!

Somehow I don't have permission to attach files... The code is pasted here:

#! /usr/bin/env python

from nutils import *

def main( nelems=32, theta_v=numpy.pi/200*numpy.array([1,16,22,25]), maxrefine=0 ):

  for theta_n in numpy.arange(theta_v.shape[0]):
    log.user('--------------------')
    log.user('domain',str(int(theta_n)),'of',theta_v.shape[0])
    theta = theta_v[theta_n]

    verts = numpy.linspace( -1, 1, 2*nelems+1 )
    domain, geom = mesh.rectilinear( [verts,verts] )
    Q = function.asarray([ [  numpy.cos(theta), -numpy.sin(theta)],
                           [  numpy.sin(theta),  numpy.cos(theta)] ])
    x ,y     = geom
    xi,gamma = (Q.T*geom[_,:]).sum(-1)

    domain = domain.trim( .5 - xi   , maxrefine=maxrefine, name = 'trimright' )
    domain = domain.trim( .5 + xi   , maxrefine=maxrefine, name = 'trimleft' )
    domain = domain.trim( .5 - gamma, maxrefine=maxrefine, name = 'trimtop' )
    domain = domain.trim( .5 + gamma, maxrefine=maxrefine, name = 'trimbottom' )
    domain.volume_check(geom)

    elem_dict = {}

    #Indexing boundary elements
    for belem in domain.boundary:
      etrans = belem.transform.lookup( domain.edict )
      ielem  = domain.edict[ etrans ]
      elem_dict.setdefault(ielem,[]).append( belem )

    #Loop over boundary elements
    for ielem, belems in log.iter( 'boundary element', elem_dict.items() ):
      velem = domain.elements[ielem]

      #Create local topology
      vdomain = topology.Topology( [velem] )
      bdomain = topology.Topology( belems )

      #Calculate centre of mass
      length = bdomain.integrate( 1., geometry=geom, ischeme='gauss1' )
      volume, centre_mass_vol = vdomain.integrate( [ 1., geom ], geometry=geom, ischeme='gauss1' )
      centre_mass = centre_mass_vol/volume

      if volume < 1e-6:
        log.user('----------')
        log.user('element',ielem)
        log.user('length',length)
        log.user('volume',volume)
        log.user('centre_mass',centre_mass)
        vol_coords, vals = vdomain.simplex.elem_eval( [geom, 1.], ischeme='bezier5', separate=True )
        bound_coords     = bdomain.simplex.elem_eval(  geom     , ischeme='bezier5', separate=True )
        with plot.PyPlot( 'theta_'+str(int(theta_n))+'_element_'+str(int(ielem)), ndigits=0 ) as plt:
          plt.mesh( vol_coords, vals, edgewidth=.1, triangulate = 'bezier')
          for line_n in numpy.arange(len(bound_coords)):
            plt.plot(bound_coords[line_n][:,0],bound_coords[line_n][:,1],c='k',linewidth=10)
          plt.axes().set_aspect('auto')
        log.user('----------')
    log.user('--------------------')

  return

util.run( main )

Documentation

Missing: an introductory document that explains fundamental concepts like topologies and functions.

add function decorator with serialization cache

When calling a decorated function all array arguments are replaced by placeholders. The function is serialized once for each unique set of placeholder arguments. Decorated functions behave exactly as undecorated functions provided that the function does not rely on values outside the function namespace, i.e. returns the same value whenever called with the same set of arguments.

Difference operation on two WithChildren elements

When constructing the complementary domain after a second trimming operation, the following error is obtained:

TypeError("unsupported operand type(s) for -: 'WithChildrenReference' and 'WithChildrenReference'",)
  File "../../nutils/nutils/topology.py", line 111, in Topology.__sub__
    refs = [ elem.reference - refmap.pop(elem.transform,elem.reference.empty) for elem in self ]

Run the code below to reproduce:

from nutils import *

def main ( ne = 10,
           nd = 2 ,
           mr = 2 ,
           r1 = 0.13,
           x1 = (0.4,0.6),
           r2 = 0.16,
           x2 = (0.6,0.3) ):

  topo, geom = mesh.rectilinear( [numpy.linspace(0,1,ne+1)]*nd )

  levelset = (geom[0]-x1[0])**2+(geom[1]-x1[1])**2-r1**2
  dom1 = topo.trim( levelset, maxrefine=mr, eps=1/1024. )
  dom2 = topo - dom1

  levelset = (geom[0]-x2[0])**2+(geom[1]-x2[1])**2-r2**2
  topo = dom1
  dom1 = topo.trim( levelset, maxrefine=mr, eps=1/1024. )
  dom3 = topo - dom1

  #Plotting
  points1, ids1 = dom1.simplex.elem_eval( [geom,1], ischeme='vtk', separate=True )
  points2, ids2 = dom2.simplex.elem_eval( [geom,2], ischeme='vtk', separate=True )
  points3, ids3 = dom3.simplex.elem_eval( [geom,3], ischeme='vtk', separate=True )
  points = points1 + points2 + points3
  ids    = ids1    + ids2    + ids3
  with plot.VTKFile('trimmed',ndigits=0) as vtkfile:
    vtkfile.unstructuredgrid( points )
    vtkfile.pointdataarray( 'ids', ids )

util.run( main )

gmesh reader

On the master branch, somewhere between release and tip, something seems to upset the gmesh reader (not fixed by some of the gmesh-fix commits). A sample code demonstrating this can be found here. Anyone else using the gmesh reader and running into similar problems?

Roadmap nutils

It would be preferable to have a roadmap for future nutils developments, including an indication of when to expect certain developments. This roadmap should not be there to nag the BDFL, but it should provide a basis for reflecting on the efforts we put in this project.

add Trilinos matrix backend

outdir property in nutilsrc stopped working

(9e3ebd2 or earlier) The output is written to public_html, despite the provision of a manually set "outdir" property in the nutilsrc file.

implement tensorial topologies

Definition of function.jump counter-intuitive

The jump as defined in function.py seems somewhat counter-intuitive. It would be preferable to define it as the opposite minus self.

Collaborative logo effort

With the nutils project online and under active development, the one thing most glaringly missing at this point is a project logo.

That I was not the only one to notice this omission was clear from a number of unsolicited suggestions, all highly appreciated, that were sent to me in the weeks past. Clearly the nutils theme inspires the designer in us. To leave everybody in on the fun it seems only fair to make this an open call for logos. I would call it a contest, but contests call for a jury and other impracticalities. Besides, as with science, I believe that openness and discussion lead to better results than blind competition.

So here we have it: I call on everybody so inclined to reply to this thread with logo ideas attached. You will need to use the github web interface to attach images. Don't worry about image quality or file formats: at this stage it is the idea that counts. If you lack the practical skills, but think you have a good idea to share then by all means do. Or if you see something with potential and have ideas for improvement, don't hesitate to build on other people's work. Be assured that by your shear participation in this list you will gain acknowledgement for the end result.

Some considerations. The project is called nutils short for Numerical Utilities. The alternative pronunciation "noodles" lends itself well for a logo, and I personally welcome artwork inspired around that theme. However, a really strong logo would combine the gimmick with a reference to the actual meaning. This is hard, and I don't have a ready idea of how to do that. But together we may think of something.

Looking forward to your ideas!

Unit test fails: test_function.TestInv2x2.test_doublegradient

I got an error while running make test

test_topology.TestConnectivityStructuredPeriodic.test_1DConnectivity ... ok
test_topology.TestConnectivityStructuredPeriodic.test_2DConnectivity ... ok
test_topology.TestConnectivityStructuredPeriodic.test_3DConnectivity ... ok
test_topology.TestStructure2D.testBoundaries ... ok
test_topology.TestStructure2D.testMesh ... ok

FAIL: test_function.TestInv2x2.test_doublegradient

Traceback (most recent call last):
File "/Users/Nitish/anaconda/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
self.test(*self.arg)
File "/Users/Nitish/nutils/tests/test_function.py", line 104, in test_doublegradient
numpy.testing.assert_array_almost_equal( fddgrad, G.eval(self.elem,self.points), decimal=2 )
File "/Users/Nitish/anaconda/lib/python2.7/site-packages/numpy/testing/utils.py", line 811, in assert_array_almost_equal
header=('Arrays are not almost equal to %d decimals' % decimal))
File "/Users/Nitish/anaconda/lib/python2.7/site-packages/numpy/testing/utils.py", line 644, in assert_array_compare
raise AssertionError(msg)
AssertionError:
Arrays are not almost equal to 2 decimals

(mismatch 1.04166666667%)
x: array([[[[[[ -2.23044871e+01, -4.72058304e+01],
[ -4.72058304e+01, -7.73985409e+01]],
...
y: array([[[[[[ -2.23045077e+01, -4.72058552e+01],
[ -4.72058552e+01, -7.73985107e+01]],
...

---

Ran 376 tests in 79.084s

FAILED (failures=1)
make: *** [test_nose] Error 1

Names/conventions in nutils

Numbering in ProductElement:

elem1, elem2 = elem0, elem1 # consistent changing throughout nutils and tests.

Do not use shared zeros when nprocs=1

add derivatives to decorated functions

distributed function evaluations

Documentation

function.py:20
Two lines seem to have disappeared here, noticed because the grammar is now wrong.
Secondly, lever=level? Good idea to pull the documentation through a spell check before the next release maybe?
ps. Notwithstanding, I very much like the documentation, and am surprised it took me a year to stumble across this issue ;-)