Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
Sfepy code to test displacement in 3D
import numpy as np
from sfepy.base.goptions import goptions
goptions['verbose'] = True
from sfepy.discrete.fem import Field
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 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
shape = (3, 3, 3)
dx = 1.
displacement = 0.05
center = np.zeros_like(shape)
mesh = gen_block_mesh(shape, np.array(shape) + 1, center,
domain = Domain('domain', mesh)
region_all = domain.create_region('region_all', 'all')
field = Field.from_args('fu', np.float64, 'vector', region_all,
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 ='dw_lin_elastic_iso(m.lam,, 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',
print region_fix_points.vertices
fix_points_BC = EssentialBC('fix_points_BC',
{'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',
print region_shift_points.vertices
shift_points_BC = EssentialBC('shift_points_BC',
{'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)
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)

This comment has been minimized.

Copy link
Owner Author

wd15 commented Sep 4, 2014

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.