Two Problems: Define different node forces and calculating cauchy stress in a post process hook function
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.
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.
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
Am Freitag, 3. März 2017 21:37:06 UTC+1 schrieb Robert Cimrman:
Hi,
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
On 03/02/2017 05:47 PM, kathrinsbriefkasten via sfepy-devel wrote: 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.
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
coorcorresponds 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:
- define a variable of the 'parameter field' kind, with the linear Lagrange basis.
- set its DOFs to your nodal values.
- 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_materi...
Regard, Kathrin
Am Freitag, 3. März 2017 21:37:06 UTC+1 schrieb Robert Cimrman:
Hi,
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
On 03/02/2017 05:47 PM, kathrinsbriefkasten via sfepy-devel wrote: 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.
participants (2)
-
kathrinsb...＠googlemail.com -
Robert Cimrman