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

> # 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.
>
>