linear_elasticity/multi_node_lcbcs.py¶

Description

Linear elasticity with multi node linear combination constraints.

Find such that:

where

The hanging nodes in 'Dependent_nodes' region are bounded to the nodes in 'Independent_nodes' region using the lcbcs (multi_node_combination) conditions

View the results using:

sfepy-view hanging_nodes.vtk -f u:wu:e -2


source code

# -*- coding: utf-8 -*-
r"""
Linear elasticity with multi node linear combination constraints.

Find :math:\ul{u} such that:

.. math::
\int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u})
= 0

where

.. math::
D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) +
\lambda \ \delta_{ij} \delta_{kl}
\;.

The hanging nodes in 'Dependent_nodes' region are bounded
to the nodes in 'Independent_nodes' region using the lcbcs
(multi_node_combination) conditions

View the results using::

sfepy-view hanging_nodes.vtk -f u:wu:e -2
"""
import numpy as nm
from sfepy import data_dir
from sfepy.base.base import output
from sfepy.mechanics.matcoefs import stiffness_from_lame

def node_map_fun(coors0, coors1):
"""
Map hanging nodes (coors0) and coincident edges (edge nodes in coors1)
on a line.
"""

ldir = coors1[-1] - coors1[0]
ldir = ldir / nm.linalg.norm(ldir)
if ldir[0] > 1e-9:
par0 = coors0[:, 0] / ldir[0]
par1 = coors1[:, 0] / ldir[0]
else:
par0 = coors0[:, 1] / ldir[1]
par1 = coors1[:, 1] / ldir[1]

idxs1 = nm.argsort(par1)
nmap = []
for par in par0:
edge = nm.where(nm.logical_and(par1[:-1] < par, par < par1[1:]))[0][0]
de = par1[edge + 1] - par1[edge]
nmap.append((idxs1[edge], idxs1[edge + 1],
1 - (par - par1[edge]) / de,
1 - (par1[edge + 1] - par) / de))

return (nm.arange(par0.shape[0]),
nm.array([k[:2] for k in nmap]),
nm.array([k[2:] for k in nmap]))

def check_lcs(out, pb, state, extend=None):
displ = state.get_state_parts()['u'].reshape((-1, 2))
coors = pb.domain.mesh.coors

hnodes_map = {12: [2, 6], 14: [2, 6], 16: [6, 10], 18: [6, 10]}

ok = True
for nd, ends in hnodes_map.items():
dy = coors[ends[1], 1] - coors[ends[0], 1]
t = (coors[nd, 1] - coors[ends[0], 1]) / dy
u_exp = displ[ends[0]] * (1 - t) + displ[ends[1]] * t
ok = ok and (nm.linalg.norm(displ[nd] - u_exp) < 1e-9)

result = {True: 'passed', False: 'failed'}[ok]
output(f'multi node linear condition test {result}!')

# tests_lcbcs hook:
nls_status = pb.get_solver().status.nls_status
if hasattr(nls_status, 'conditions'):
nls_status.conditions.append(int(not ok))

return out

filename_mesh = data_dir + '/meshes/2d/special/hanging_nodes.mesh'

options = {
'post_process_hook': check_lcs,
}

regions = {
'Omega': 'all',
'Omega1': 'cells of group 1',
'Omega2': 'cells of group 2',
'Inside_nodes': ('vertices in (x > 0.199) & (x < 0.201)', 'vertex'),
'Independent_nodes': ('r.Inside_nodes *v r.Omega1', 'vertex'),
'Dependent_nodes': ('r.Inside_nodes -v r.Independent_nodes', 'vertex'),
'Left': ('vertices in (x < 0.001)', 'facet'),
'Right': ('vertices in (x > 0.299)', 'facet'),
}

materials = {
'solid': ({'D': stiffness_from_lame(dim=2, lam=1e1, mu=1e0)},),
}

fields = {
'displacement': ('real', 'vector', 'Omega', 1),
}

integrals = {
'i': 1,
}

variables = {
'u': ('unknown field', 'displacement'),
'v': ('test field', 'displacement', 'u'),
}

ebcs = {
'Fixed': ('Left', {'u.all': 0.0}),
'Displaced': ('Right', {'u.0': 0.05, 'u.1': -0.05}),
}

functions = {
'node_map_fun': (node_map_fun,),
}

lcbcs = {
'hanging': (['Dependent_nodes', 'Independent_nodes'], {'u.all': 'u.all'},
'node_map_fun', 'multi_node_combination', None),
}

equations = {
'balance_of_forces' :
"""dw_lin_elastic.i.Omega(solid.D, v, u) = 0""",
}

solvers = {
'ls': ('ls.auto_direct', {}),
'newton': ('nls.newton', {
'i_max': 1,
'eps_a': 1e-10,
}),
}