linear_elasticity/mixed_mesh.py

Description

A linear elastic beam loaded with a continuous force. The FE mesh consists of hexehedral and tetrahedral elements.

The displacement at the beam end is compared to the reference solution calculated on the homogeneous hexahedral mesh.

Running the simulation:

sfepy-run sfepy/examples/linear_elasticity/mixed_mesh.py

Viewing the results:

sfepy-view beam_h* -f u:s0:wu:e:p0 u:s1:wu:e:p0 --camera-position="1.2,-0.6,0.1,0.4,0.1,-0.1,-0.2,0.1,1"

source code

r"""
A linear elastic beam loaded with a continuous force. The FE mesh consists
of hexehedral and tetrahedral elements.

The displacement at the beam end is compared to the reference
solution calculated on the homogeneous hexahedral mesh.

Running the simulation::

    sfepy-run sfepy/examples/linear_elasticity/mixed_mesh.py

Viewing the results::

    sfepy-view beam_h* -f u:s0:wu:e:p0 u:s1:wu:e:p0 --camera-position="1.2,-0.6,0.1,0.4,0.1,-0.1,-0.2,0.1,1"
"""
from __future__ import absolute_import
from sfepy import data_dir
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
from sfepy.base.base import Struct
import numpy as nm


def reference_solution(pb):
    from sfepy.discrete import Problem

    conf = pb.conf.copy()
    conf.filename_mesh = data_dir + '/meshes/3d/beam_h7.mesh'
    for k in list(conf.regions.keys()):
        if '_t' in k:
            del conf.regions[k]

    for k in list(conf.fields.keys()):
        if '_t' in k:
            del conf.fields[k]

    for k in list(conf.variables.keys()):
        if '_t' in k:
            del conf.variables[k]

    conf.equations = {
        'balance_of_forces':
            """dw_lin_elastic.i.Omega(solid.D, v, u)
             = dw_surface_ltr.i.Top(force.val, v)"""
    }

    rpb = Problem.from_conf(conf)
    rpb.set_output_dir(pb.output_dir)

    return rpb.solve().get_state_parts()['u']


def get_force(ts, coors, mode=None, **kwargs):
    if mode == 'qp':
        F = 1e3

        val = nm.zeros_like(coors)[..., None]
        val[:, 2, 0] = -coors[:, 0] / 0.7 * F

        return {'val': val}


def post_proces(out, pb, state, extend=False):
    c1 = pb.domain.regions['Omega'].cells
    c2 = pb.domain.regions['Omega_t'].cells
    ed = pb.domain.regions['Edge'].vertices

    S = nm.zeros((pb.domain.cmesh.n_el, 1, 6, 1), dtype=nm.float64)

    S[c1] = pb.evaluate('ev_cauchy_stress.i.Omega(solid.D, u)',
                        mode='el_avg')
    S[c2] = pb.evaluate('ev_cauchy_stress.i.Omega_t(solid.D, u_t)',
                        mode='el_avg')

    out['stress'] = Struct(name='out', data=S, mode='cell')

    u_nd = out['u'].data[ed[0], :]
    ref_u_nd = reference_solution(pb).reshape(-1, 3)[ed[0], :]

    du = nm.linalg.norm((ref_u_nd - u_nd) / nm.linalg.norm(ref_u_nd))
    print(f'Relative difference with respect to the reference solution: {du}')

    return out


filename_mesh = data_dir + '/meshes/3d/beam_h5t12.mesh'

options = {
    'post_process_hook': 'post_proces',
}

regions = {
    'Omega_': 'all',
    'Omega': ('cells of group 1', 'cell', None, {'vertices_from': 'Omega_'}),
    'Omega_t': ('cells of group 2', 'cell', None, {'vertices_from': 'Omega_'}),
    'Left': ('vertices in (x < 0.01)', 'facet'),
    'Right': ('vertices in (x > 0.69)', 'facet'),
    'Top_': ('vertices in (z > 0.09)', 'facet'),
    'Bottom': ('vertices in (z < 0.01)', 'facet'),
    'Edge': ('r.Bottom *v r.Right', 'vertex'),
    'Top': ('r.Top_ *f r.Omega', 'facet'),
    'Top_t': ('r.Top_ *f r.Omega_t', 'facet'),
}

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

materials = {
    'solid': ({'D': stiffness_from_youngpoisson(dim=3, young=1e6, poisson=0.3)},),
    'force': 'get_force',
    'vforce': ({'.val' : [0.0, 0.0, -7]},),
}

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

integrals = {
    'i': 2,
}

variables = {
    'u': ('unknown field', 'displacement', 0),
    'v': ('test field', 'displacement', 'u'),
    'u_t': ('unknown field', 'displacement_t', 0),
    'v_t': ('test field', 'displacement_t', 'u_t'),
}

ebcs = {
    'Fixed': ('Left', {'u.all' : 0.0}),
}

equations = {
    'balance_of_forces':
    """dw_lin_elastic.i.Omega(solid.D, v, u)
     + dw_lin_elastic.i.Omega_t(solid.D, v_t, u_t)
     = dw_surface_ltr.i.Top(force.val, v)
     + dw_surface_ltr.i.Top_t(force.val, v_t)"""
}

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