Sfepy code to test periodic displacements

sfepy_lcbcs.py

import numpy as np
from sfepy.base.goptions import goptions
goptions['verbose'] = False
from sfepy.discrete.fem import Field
try:
    from sfepy.discrete.fem import FEDomain as Domain
except ImportError:
    from sfepy.discrete.fem import Domain
from sfepy.discrete import (FieldVariable, Material, Integral, Function,
                            Equation, Equations, Problem)
from sfepy.terms import Term
from sfepy.discrete.conditions import Conditions, EssentialBC, PeriodicBC
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
import sfepy.discrete.fem.periodic as per
from sfepy.discrete import Functions
from sfepy.mesh.mesh_generators import gen_block_mesh
from sfepy.base.base import output
from sfepy.discrete.conditions import LinearCombinationBC
from sfepy.postprocess.viewer import Viewer


output.set_output(quiet=True)

shape = (3, 3)
dx = 1.

center = np.zeros_like(shape)
mesh = gen_block_mesh(shape, np.array(shape) + 1, center,
                      verbose=False)

domain = Domain('domain', mesh)

region_all = domain.create_region('region_all', 'all')

field = Field.from_args('fu', np.float64, 'vector', region_all,
                        approx_order=2)

u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')

m = Material('m', lam=0.576, mu=0.256)
integral = Integral('i', order=3)

t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u)',
              integral, region_all, m=m, v=v, u=u)
eq = Equation('balance_of_forces', t1)
eqs = Equations([eq])

## periodic boundary condition for y-plane

eps = 1e-3 * dx
min_xy, max_xy = domain.get_mesh_bounding_box()

def cutplane(coords, X, index):
    c = coords[:, index]
    return (c > (X - eps)) & (c  < (X + eps))

yplus_ = lambda c, domain=None: np.where(cutplane(c, max_xy[1], 1))[0]
yminus_ = lambda c, domain=None: np.where(cutplane(c, min_xy[1], 1))[0]

yplus = Function('yplus', yplus_)
yminus = Function('yminus', yminus_)

region_top = domain.create_region('region_y_plus',
                                  'vertices by yplus',
                                  'facet',
                                  functions=Functions([yplus]))


region_bottom = domain.create_region('region_y_minus',
                                     'vertices by yminus',
                                     'facet',
                                     functions=Functions([yminus]))

match_x_line = Function('match_x_line', per.match_x_line)

epbc = PeriodicBC('periodic_y', 
                  [region_top, region_bottom],
                  {'u.all': 'u.all'},
                  match='match_x_line')

epbcs, functions = Conditions([epbc]), Functions([match_x_line])

## fixed displacement boundary conditions

xminus_ = lambda c, domain=None: np.where(cutplane(c, min_xy[0], 0))[0]
xminus = Function('xminus', xminus_)


def fix_x_points_(c, domain=None):
    mask = cutplane(c, min_xy[0], 0) & (cutplane(c, min_xy[1], 1) | cutplane(c, max_xy[1], 1))
    return np.where(mask)[0]

fix_x_points = Function('fix_x_points', fix_x_points_)
region_fix_points = domain.create_region('region_fix_points',
                                         'vertices by fix_x_points',
                                         'vertex',
                                         functions=Functions([fix_x_points]))

fix_points_BC = EssentialBC('fix_points_BC',
                            region_fix_points,
                            {'u.0' : 0.0, 'u.1': 0.0})
ebcs = Conditions([fix_points_BC])

## displaced boundary conditions

xplus_ = lambda c, domain=None: np.where(cutplane(c, max_xy[0], 0))[0]
xplus = Function('xplus', xplus_)

def shift_x_points_(c, domain=None):
    mask = cutplane(c, max_xy[0], 0) & (cutplane(c, min_xy[1], 1) | cutplane(c, max_xy[1], 1))
    return np.where(mask)[0]

shift_x_points = Function('shift_x_points', shift_x_points_)
region_shift_points = domain.create_region('region_shift_points',
                                           'vertices by shift_x_points',
                                           'vertex',
                                           functions=Functions([shift_x_points]))

shift_points_BC = EssentialBC('shift_points_BC',
                              region_shift_points,
                              {'u.0' : 0.3, 'u.1': 0.0})

ebcs = Conditions([fix_points_BC, shift_points_BC])

## linear combinations boundary conditions 

xplus_ = lambda c, domain=None: np.where(cutplane(c, max_xy[0], 0))[0]
xplus = Function('xplus', xplus_)


region_x_plus = domain.create_region('region_x_plus',
                                     'vertices by xplus',
                                     'facet',
                                     functions=Functions([xplus]))
region_x_minus = domain.create_region('region_x_minus',
                                      'vertices by xminus',
                                      'facet',
                                      functions=Functions([xminus]))
match_y_line = Function('match_y_line', per.match_y_line)
def shift_(ts, coors, region):
    return 0.3 * np.ones_like(coors[:, 0])
shift = Function('shift', shift_)
lcbc = LinearCombinationBC('lcbc',
                           [region_x_plus, region_x_minus],
                           {'u.0' : 'u.0'},
                           match_y_line,
                           'shifted_periodic',
                           arguments=(shift,))
        
lcbcs = Conditions([lcbc])

## solve equations

ls = ScipyDirect({})

pb = Problem('elasticity', equations=eqs, auto_solvers=None)
pb.save_regions_as_groups('regions')

pb.time_update(ebcs=ebcs, epbcs=epbcs, lcbcs=lcbcs, functions=functions)

ev = pb.get_evaluator()
nls = Newton({}, lin_solver=ls,
             fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix)

pb.set_solvers_instances(ls, nls)

vec = pb.solve()

u = vec.create_output_dict()['u']
print u.data[:4,0]
print u.data[-4:,0]

print u.data[:,0]


pb.save_state('dump.vtk', vec)
view = Viewer('dump.vtk')
view(vector_mode='warp_norm', rel_scaling=2, is_scalar_bar=True, is_wireframe=True)

The above code tests periodic displacements in Sfepy. The code was run with https://github.com/rc/sfepy/tree/a9fdd2b4d915d2d07407d5ba246734700e0c5094, but with the following modification:

diff --git a/sfepy/discrete/fem/lcbc_operators.py b/sfepy/discrete/fem/lcbc_operators.py
index af8bf6a..dd5be0d 100644
--- a/sfepy/discrete/fem/lcbc_operators.py
+++ b/sfepy/discrete/fem/lcbc_operators.py
@@ -358,7 +358,7 @@ class ShiftedPeriodicOperator(LCBCOperator):

         num = len(ia)

-        ones = nm.ones(num, dtype=nm.float64)
+        ones = -nm.ones(num, dtype=nm.float64)
         n_dofs = [variables.adi.n_dof[name] for name in self.var_names]
         mtx = sp.coo_matrix((ones, (meq, seq)), shape=n_dofs)

The modification seems to make the results better in terms of displacement direction, but not correct.

The mesh is just a 2D square with 3x3 elements (each element has side of length 1). The boundary conditions are

• fix the left most top and bottom node in both the x and y directions
• the top and bottom plane have regular periodic displacements in both the x and y-directions
• the left and right planes have a periodic offset of 0.3 in the x-direction

The periodic offset is set with line19, the shift_ function. I am assuming that the way the shift_ function is written is correct.

At the end of the example above, the displacements in the x and y-directions are displayed using print u.data[:4,0] and the right side print u.data[-4:,0]. The results are

[ 0.         0.1212875  0.1212875  0.       ]

for the left side and

[ 0.15564665  0.1787125   0.1787125   0.15564665]

for the right side. The top and bottom left side displacements are 0, as expected since the points are fixed. The top and bottom right side points have the same displacement (0.1556), but that isn't 0.3, which is what I would expect since that is what I'm assuming the shift_ function is doing. Moreover, the two center nodes on each side have an offset of only 0.178 - 0.121 = 0.057 (not 0.3). Furthermore, the offset is different between the center nodes and the corner nodes, this should at least be the same down the entire plane.

How do I make this work? How do I get the offset to be the same for each left and right node and how do I get that value to be 0.3.

Note that I haven't fixed the periodic displacement between the left and right side in the y-direction. The offsets in the y-direction are quite small so I haven't worried about this issue yet.

I have made some fixes:

• use per.match_x_line instead of per.match_y_plane
• use per.match_y_line instead of per.match_x_plane - the "plane" functions are for 3D.

I changed the approximation order to one, to have only vertex displacements for simplicity of debugging and smaller vectors.
With that, and unmodified sfepy, I get:

[ 0.         -0.16176654 -0.16176654  0.        ]
[-0.02353307  0.13823346  0.13823346 -0.02353307]

and

In [11]: 0.13823346 - -0.16176654
Out[11]: 0.3

The first and last values are given by other BCs (fixed left corner nodes prevent using LCBC in corners). So this looks ok to me(?)

With approx_order two:

[ 0.        -0.1170099 -0.1170099  0.       ]
[ 0.06598021  0.1829901   0.1829901   0.06598021]

and

In [12]: 0.1829901 - -0.1170099
Out[12]: 0.3

Note that u = vec.create_output_dict()['u'] by default strips (throws away) the higher order DOFs.

If you want to have the shift 0.3 also in the corner nodes, just fix by EBC the u.0 displacement in the right-hand side corners to 0.3. Then I get:

[  0.00000000e+00   3.08781635e-16   2.21055642e-16   0.00000000e+00]
[ 0.3  0.3  0.3  0.3]

Btw. if I then remove the LCBC:

[ 0.          0.15566122  0.15566122  0.        ]
[ 0.3         0.14433878  0.14433878  0.3       ]

So above the all 0.3 were really enforced by the LCBC.

When I change the shape to e.g. (100, 100), the corner node influence diminishes:

[ 0.         -0.00697834 -0.00697834 -0.00697834]
[ 0.29302166  0.29302166  0.29302166  0.28604332]

=> It looks ok - do you agree?