Sfepy code to test displacement in 3D

gistfile1.py

import numpy as np
from sfepy.base.goptions import goptions
goptions['verbose'] = True
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
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.discrete import Functions
from sfepy.mesh.mesh_generators import gen_block_mesh
from sfepy.base.base import output
from sfepy.postprocess.viewer import Viewer


output.set_output(quiet=False)

shape = (3, 3, 3)
dx = 1.
displacement = 0.05

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

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])

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))

## 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)) & \
    (cutplane(c, min_xy[2], 2) | cutplane(c, max_xy[2], 2))
    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]))
print region_fix_points.vertices
fix_points_BC = EssentialBC('fix_points_BC',
                            region_fix_points,
                            {'u.0' : 0.0, 'u.1': 0.0, 'u.2': 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)) & \
    (cutplane(c, min_xy[2], 2) | cutplane(c, max_xy[2], 2))
    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]))
print region_shift_points.vertices
shift_points_BC = EssentialBC('shift_points_BC',
                              region_shift_points,
                              {'u.0' : displacement, 'u.1': 0.0, 'u.2': 0.0})

ebcs = Conditions([fix_points_BC, shift_points_BC])


## solve equations

ls = ScipyDirect({})

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

pb.time_update(ebcs=ebcs)

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']

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 code above runs an Sfepy FE simulation in 3D. The boundary conditions are

• fix the corner right side nodes in all directions
• displace the corner left side nodes in the x-direction and fix in y and z.

The result from this calculation looks very distorted and not symmetric at all.

What is going wrong?