Two Problems: Define different node forces and calculating cauchy stress in a post process hook function

[kathrinsb...@googlemail.com]

Hi,

I am running the following FE simulation example (linear_elastic_tractions) as a python program.

The first question is how to change the definition of the pressure traction at the surface. In my problem I have different force values for each node in the top/bottom region?

The second question is about how to calculate the stress results. I extendes the traction example with a post_process_hook. My program is as follows:

# coding: utf-8
from __future__ import print_function
from __future__ import absolute_import
from argparse import ArgumentParser
import numpy as nm
import time

import sys
sys.path.append('.')

from sfepy.base.base import IndexedStruct
from sfepy.discrete import (FieldVariable, Material, Integral,
                            Equation, Equations, Problem, Function)
from sfepy.discrete.fem import Mesh, FEDomain, Field
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.postprocess.viewer import Viewer
from sfepy.mechanics.matcoefs import stiffness_from_lame

def linear_tension(ts, coor, mode=None, **kwargs):
    if mode == 'qp':
        val = nm.tile(1.0, (coor.shape[0], 1, 1))
        return {'val' : val}

def post_process(out, problem, state, extend=False):
    from sfepy.base.base import Struct
    # Cauchy strain averaged in elements.
    strain = problem.evaluate('ev_cauchy_strain.2.Omega(u)',
                              mode='el_avg')
    out['cauchy_strain'] = Struct(name='output_data',
                                  mode='cell', data=strain,
                                  dofs=None)
    # Cauchy stress averaged in elements.
    stress = problem.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)',
                              mode='el_avg')
    out['cauchy_stress'] = Struct(name='output_data',
                                  mode='cell', data=stress,
                                  dofs=None)
    return out

def cube_displacement_medium():
    from sfepy import data_dir

    parser = ArgumentParser()
    parser.add_argument('--version', action='version', version='%(prog)s')
    parser.add_argument('-s', '--show',
                      action="store_true", dest='show',
                      default=True, help=help['show'])
    options = parser.parse_args()

    mesh = Mesh.from_file(data_dir + '/meshes/3d/cube_medium_tetra.mesh')
    print("Mesh number of nodes: {}, number of tetrahedrons: {}".format(mesh.n_nod, mesh.n_el))

    domain = FEDomain('domain', mesh)
    domain.refine()

    omega = domain.create_region('Omega', 'all') 
    left = domain.create_region('Left',
                                  'vertices in (x < -0.1)',
                                  'facet')
    right = domain.create_region('Right',
                               'vertices in (x > 0.3)',
                               'facet')
    field = Field.from_args('displacement', nm.float64, 3, omega, approx_order=1)
    u = FieldVariable('u', 'unknown', field, order=0)
    v = FieldVariable('v', 'test', field, primary_var_name='u')
   
    solid = Material('solid', values={'D': stiffness_from_lame(dim=3, lam=5.769, mu=3.846)})
    load = Material('load', function=linear_tension)
   
    integral = Integral('i', order=2)
   
    t1 = Term.new('dw_lin_elastic(solid.D, v, u)', integral, omega, solid=solid, v=v, u=u)
    t2 = Term.new('dw_surface_ltr(load.val, v)', integral, right, load=load, v=v)
    eq = Equation('elasticity', t1 - t2)
    eqs = Equations([eq])
    fixb = EssentialBC('fixb', left, {'u.all' : 0.0}) 
    fixt = EssentialBC('fixt', right, {'u.[1,2]' : 0.0})
    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({ 'i_max'      : 1,
                   'eps_a'      : 1e-10,
                   'eps_r'      : 1.0,
                   'macheps'   : 1e-16,
                   # Linear system error < (eps_a * lin_red).
                   'lin_red'    : 1e-2,
                   'ls_red'     : 0.1,
                   'ls_red_warp' : 0.001,
                   'ls_on'      : 1.1,
                   'ls_min'     : 1e-5,
                   'check'     : 0,
                   'delta'     : 1e-6, }, lin_solver=ls, status=nls_status)

    pb = Problem('linear_elastic_traction', equations=eqs, nls=nls, ls=ls)
    pb.save_regions_as_groups('regions')

    pb.time_update(ebcs=Conditions([fixb, fixt]))
    vec = pb.solve()
    print(nls_status)

    # calculate stress
    pb.save_state('linear_elastic_traction_cylinder.vtk', vec, post_process_hook=post_process)

    if options.show:
        view = Viewer('linear_elastic_traction_cylinder.vtk')
        view(vector_mode='warp_norm', rel_scaling=2,
            is_scalar_bar=True, is_wireframe=True)

if __name__ == '__main__':
    cube_displacement_medium()


When I am running the program I get the following error from the post process function:


sfepy: ev_cauchy_stress.2.Omega(solid.D, u)
Traceback (most recent call last):
  File "/home/experimental/fem/linear_elastic_traction_tetraeder.py", line 440, in <module>
    runtime_small_cube, nodes_small_cube = cube_displacement_medium()
  File "/home/experimental/fem/linear_elastic_traction_tetraeder.py", line 345, in cube_displacement_medium
    pb.save_state('linear_elastic_traction_cylinder.vtk', vec, post_process_hook=post_process)
  File "/bin/python27/lib/python2.7/site-packages/sfepy-2016.3_git_f99ff34eeb05ea11e61aa3639c47537daba7cc64-py2.7-linux-x86_64.egg/sfepy/discrete/problem.py", line 773, in save_state
    out = post_process_hook(out, self, state, extend=extend)
  File "/home/experimental/fem/linear_elastic_traction_tetraeder.py", line 85, in post_process
    mode='el_avg')
  File "/bin/python27/lib/python2.7/site-packages/sfepy-2016.3_git_f99ff34eeb05ea11e61aa3639c47537daba7cc64-py2.7-linux-x86_64.egg/sfepy/discrete/problem.py", line 1262, in evaluate
    verbose=verbose, **kwargs)
  File "/bin/python27/lib/python2.7/site-packages/sfepy-2016.3_git_f99ff34eeb05ea11e61aa3639c47537daba7cc64-py2.7-linux-x86_64.egg/sfepy/discrete/problem.py", line 1203, in create_evaluable
    kwargs=kwargs)
  File "/bin/python27/lib/python2.7/site-packages/sfepy-2016.3_git_f99ff34eeb05ea11e61aa3639c47537daba7cc64-py2.7-linux-x86_64.egg/sfepy/discrete/evaluate.py", line 213, in create_evaluable
    user=extra_args, verbose=verbose)
  File "/bin//python27/lib/python2.7/site-packages/sfepy-2016.3_git_f99ff34eeb05ea11e61aa3639c47537daba7cc64-py2.7-linux-x86_64.egg/sfepy/discrete/equations.py", line 66, in from_conf
    materials, integrals, user=user)
  File "/bin/python27/lib/python2.7/site-packages/sfepy-2016.3_git_f99ff34eeb05ea11e61aa3639c47537daba7cc64-py2.7-linux-x86_64.egg/sfepy/discrete/equations.py", line 751, in from_desc
    terms.assign_args(variables, materials, user)
  File "/bin/python27/lib/python2.7/site-packages/sfepy-2016.3_git_f99ff34eeb05ea11e61aa3639c47537daba7cc64-py2.7-linux-x86_64.egg/sfepy/terms/terms.py", line 295, in assign_args
    term.assign_args(variables, materials, user)
  File "/bin/python27/lib/python2.7/site-packages/sfepy-2016.3_git_f99ff34eeb05ea11e61aa3639c47537daba7cc64-py2.7-linux-x86_64.egg/sfepy/terms/terms.py", line 506, in assign_args
    % arg_name)
ValueError: material argument solid not found!

It seems that the defined material is not known in the function. Can anyone help me? Thank you.

[Robert Cimrman]

Hi,

On 03/02/2017 05:47 PM, kathrinsbriefkasten via sfepy-devel wrote:
> Hi,
>
> I am running the following FE simulation example (linear_elastic_tractions)
> as a python program.
>
> The first question is how to change the definition of the pressure traction
> at the surface. In my problem I have different force values for each node
> in the top/bottom region?

You can define the traction by a function, as in the example you use. Just
provide a different value for each coordinate in `coor`.

> The second question is about how to calculate the stress results. I
> extendes the traction example with a post_process_hook. My program is as
> follows:

First, you do not need to use the post-process hook in script-like examples.
See for example [1] for how to compute stress and add it to output.

To fix the error, try passing copy_materials=False to problem.evaluate().

I am travelling and will be without a computer for the next week, so I cannot
help you more for now.

r.

[1] http://sfepy.org/doc-devel/examples/linear_elasticity/modal_analysis.html

[kathrinsb...@googlemail.com]

Hi Robert,

thank you for the answer. 

The definition of different values for the traction in a function works. As I found in the documentation the length of `coor` corresponds to the quadrature points. Did I understand that correctly?

My problem is that I have values for each facet (triangle) or node but not for the quadrature points. Is there a function to map values from facets/nodes to quadrature points?

Regard,

Kathrin 

[Robert Cimrman]

Hi Kathrin,

On 03/12/2017 08:04 PM, kathrinsbriefkasten via sfepy-devel wrote:
> Hi Robert,
>
> thank you for the answer.
>
> The definition of different values for the traction in a function works. As
> I found in the documentation the length of `coor` corresponds to the
> quadrature points. Did I understand that correctly?

Yes, the quadrature points in all physical elements.

> My problem is that I have values for each facet (triangle) or node but not
> for the quadrature points. Is there a function to map values from
> facets/nodes to quadrature points?

You could proceed as follows:

1. define a variable of the 'parameter field' kind, with the linear Lagrange basis.
2. set its DOFs to your nodal values.
3. proceed as in [1] - use ev_volume_integrate term to interpolate the values
into QP. Make sure to use the same integral as in your equations, to have the
same quadrature points.

Does it help, or do you need more details? In the latter case, send me your
example script, so that I could run it.

r.

[2]
http://sfepy.org/doc-devel/examples/diffusion/poisson_field_dependent_material.html