Calculate total force from cauchy stress
by Patricia Garcia Cañadilla
Dear Robert,
I have a 2D body of hyperelastic material which contracts and I would like
to compute the total force developed by the body from the cauchy stress. I
am trying to follow some of your indications I found in this group, but I
still couldn't make it works. Could you please help me to fix the problem?
I am getting the following error:
key = (region.name, integral.order, integration)
AttributeError: 'dict' object has no attribute 'name'
I am trying to do the following, inside stress_strain post-processing
function:
def stress_strain(out, problem, state, extend = False ):
from sfepy.base.base import Struct
from sfepy.mechanics.tensors import StressTransform
from sfepy.mechanics.tensors import transform_data
from sfepy.discrete.common.mappings import get_normals
ev = problem.evaluate
field = problem.fields['displacement']
region = problem.domain.regions['Gamma']
integral = problem.integrals['i2']
n = get_normals(field,integral,regions)
stress = ev('dw_tl_fib_a.1.Omega(f1.fmax, f1.eps_opt, f1.s, f1.fdir,
f1.act, v, u )',mode='qp', term_mode= 'stress');
F = ev('ev_def_grad.1.Omega(u)',mode='el_avg');
transform = StressTransform(F)
Cstress = transform.get_cauchy_from_2pk(stress)
T = Cstress*n;
Force = ev('ev_surface_integrate.2.Gamma(T)')
And here it is part of the problem configuration file.
fields = {
'displacement': ('real', 'vector', 'Omega', 1),
}
materials = {
'solid' : (None, 'get_elastic_pars'),
'load' : (None, 'linear_tension'),
'f1' : 'get_pars_fibres1',
}
variables = {
'u': ('unknown field', 'displacement', 0),
'v': ('test field', 'displacement', 'u'),
}
regions = {
'Omega' : 'all',
'Fix1' : ('vertices in x < %.10f' % (fix_point + eps2), 'facet'),
'Fix2' : ('vertices in x > %.10f' % (fix_point - eps2), 'facet'),
'Fix' : ('r.Fix1 *v r.Fix2', 'facet'),
'Gamma' : ('vertices of surface','edge'),
}
ebcs = {
'fixb' : ('Fix', {'u.all' : 0.0}),
}
integrals = {
'i1' : ('v', 1),
'i2' : ('s', 2),
}
9 months
- 4
- 5
Re: Torque
by Robert Cimrman
In examples/large_deformation/hyperelastic.py a rotation by displacements is applied. By using a similar function the vectors defining the force couples could be defined for dw_surface_ltr (IMHO). Does it make sense?
r.
----- Reply message -----
From: "Andre Smit" <freev...(a)gmail.com>
To: <sfepy...(a)googlegroups.com>
Subject: Torque
Date: Sat, Dec 18, 2010 05:10
What is the best way to apply a torque load to a model?
--
Andre
--
You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To post to this group, send email to sfepy...(a)googlegroups.com.
To unsubscribe from this group, send email to sfepy-devel...(a)googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
1 year, 7 months
- 4
- 7
Correct element orientation for tetrahedral mesh
by Kathrin Sobe
Hi Robert,
I was able to load a mesh from a nastran .bdf file. The mesh is build of tetrahedrons (nodes: 12593, tetrahedrons: 310106).
When I create the Domain it stops with the following error:
sfepy: warning: bad element orientation, trying to correct...
sfepy: warning: bad element orientation, trying to correct...
Traceback (most recent call last):
runtime, nodes = tool_displacement()
domain = FEDomain('domain', mesh)
self.fix_element_orientation()
raise RuntimeError('elements cannot be oriented! (%s)' % key)
RuntimeError: elements cannot be oriented! (3_4)
What is a correct element orientation for tetrahedral elements? The Nastran file was generated automatically.
Thanks for your help,
Kathrin
2 years, 1 month
- 2
- 1
Fwd: SfePy - Electro-thermo-mechanics problem
by Robert Cimrman
Hi Manuel,
I am forwarding your message to the mailing list, because I could not reach you
using your address (Connection timed out error). I will answer to the list as well.
r.
-------- Forwarded Message --------
Subject: SfePy - Electro-thermo-mechanics problem
Date: Mon, 28 Aug 2017 12:31:33 +0200
From: Manuel Feuchter <manuel.feuchter(a)twt-gmbh.de>
To: cimrman3(a)ntc.zcu.cz
Hey Robert,
first of all: SfePy is an amazing project! I am getting to like SfePy more and
more :)
I tried to start a thread at
https://mail.python.org/mm3/archives/list/sfepy@python.org
But unfortunately am I not able to. I recieve a "Server Error (500)" . So I try
it via Email :)
Concerning SfePy: I get stuck with a few problems/questions:
Background: I would like to set up a subsequent coupled system (so not fully
coupled):
1. electric field evolves an electric current (at the end of the day, I would
like to apply the electric field as an time (time step) dependent function.)
2. electric current results in liberated heat by joule heating
3. the liberated heat drives the thermal expansion und results in thermal stresses
# 1:
------
Concerning example:
http://sfepy.org/doc-devel/examples/multi_physics/thermal_electric.html
(multi_physics/thermal_electric.py)
Equation: dw_laplace.2.Omega( m.electric_conductivity, psi, phi ) = 0
Your comment states:
"""
First solve the stationary electric conduction problem. Then use its
results to solve the evolutionary heat conduction problem.
"""
Alright. However, just to understand it comprehensively? When you say you solve
the "stationary electric conduction problem". Does this imply actually two
equations?
$(1) \textbf{E} = -\nabla \phi$\\
$(2) \textbf{j} = \frac{1}{\varrho}\cdot\textbf{E}$
Is this right?
# 2:
------
Assuming I got it right, I would like to solve the electric conduction problem
in terms of a quasi-static circumstance (external electric potential and
internal resistivity will change).
This means I would like to solve the electric conduction problem for every time
step that I run the simulation and hence use the result (of $\textbf{j}$) for
every next time step.
I tried so, however I must have gotten sth wrong. The apparent electric field
is not used?! Am I right? And how am I able to achieve so?
(please find my code below)
# 3:
------
For the heat conduction equation: ... = dw_electric_source.2.Omega(
m.electric_conductivity, s, phi )
Here, for the source term $\textbf{p} = \varrho \cdot \textbf{j}^2$ is solved?
And how does the specific heat gets involved in your equation? Is the source
term divided by the thermal conductivity [W/mK]?
Does specific_heat to be named exactly like this to be involved in the
'dw_volume_dot.2.Omega( s, dT/dt )' term?
# 4:
------
How am I able to apply a time (time step) dependent function for the electric
potential?!
Looking forward for your reply!
Cheers and kind regards from Germany
Manuel
--------------------------------------------------------
def stress_strain_ext(out, pb, state, extend=False):
"""
Calculate and output strain and stress for given displacements.
"""
from sfepy.base.base import Struct
from sfepy.mechanics import tensors
ev = pb.evaluate
strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.2.Omega(m.mech_solid, u)', mode='el_avg')
vms = tensors.get_von_mises_stress(stress.squeeze())
vms.shape = (vms.shape[0], 1, 1, 1)
out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain, dofs=None)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress, dofs=None)
out['von_mises_stress'] = Struct(name='output_data', mode='cell',
data=vms, dofs=None)
return out
#! Presettings
#! -------
t0 = 0.0
t1 = 100.0
n_step = 11
# material parameters and constants
lam = 10.0
mu = 5.0
thermal_expandability = 1.25e-5
T0 = 20.0 # Background temperature.
specific_heat = 1.2
thermal_kappa_xyz = 2.0
electric_sigma_xyz = 1.5
#! Options
#! -------
options = {
'ts': 'ts',
'nls': 'newton',
'ls': 'ls',
'output_dir' : output_dir,
'post_process_hook' : 'stress_strain_ext',
}
#! Regions
#! -------
regions = {
'Omega': 'all',
'Omega_Cell_Surface': ('vertices of surface', 'facet'),
'Mech_Mount': ('vertices in (x < 0.05)', 'facet'),
'Mech_Front': ('vertices in (x > 0.95)', 'facet'),
'Therm_inlet': ('vertices in (x < 90.05) & (x > 85.95) & (y > 9.95)',
'facet'),
'Therm_outlet': ('vertices in (x < 0.05)', 'facet'),
'Elec_Phi_0' : ('vertices in (x < 37.05) & (x > 32.95) & (y > 9.95)',
'facet'),
'Elec_Phi_1' : ('vertices in (x < 90.05) & (x > 85.95) & (y < 0.05)',
'facet'),
}
#! Materials
#! ---------
eye_sym = np.array([[1], [1], [1], [0], [0], [0]], dtype=np.float64)
materials = {
'm': ({
'thermal_conductivity': thermal_kappa_xyz,
'electric_conductivity': electric_sigma_xyz,
'mech_solid' : stiffness_from_lame(3, lam=lam, mu=mu),
'alpha_therm' : (3.0 * lam + 2.0 * mu) * thermal_expandability * eye_sym
},),
'mech_load': ({'.val': [0.0, 0.0, 0.0]},),
}
#! Fields
#! ------
fields = {
'displacement': ('real', 'vector', 'Omega', 1),
'temperature': ('real', 1, 'Omega', 1),
'potential' : ('real', 1, 'Omega', 1),
}
#! Variables
#! ---------
variables = {
'u': ('unknown field', 'displacement', 2),
'v': ('test field', 'displacement', 'u'),
'T': ('unknown field', 'temperature', 0, 1), # 1 means, set history "on"
's': ('test field', 'temperature', 'T'),
'phi': ('unknown field', 'potential', 1, 1),
'psi': ('test field', 'potential', 'phi'),
}
#! Initial Conditions
#! -------------------
ics = {
'ic' : ('Omega', {'u.all' : 0.0,
'T.0': 0.0,
'phi.all': 0.0,
}),
}
#! Boundary Conditions
#! -------------------
ebcs = {
'u0': ('Mech_Mount', {'u.all': 0.0}),
't0': ('Mech_Mount', {'T.0': 0.0}),
't1': ('Mech_Front', {'T.0': 0.0}),
'p0' : ('Elec_Phi_0', {'phi.0' : 0.0}),
'p1' : ('Elec_Phi_1', {'phi.0' : 100.0}),
}
#! Equations
#! ---------
equations = {
'potential' :
"""dw_laplace.2.Omega( m.electric_conductivity, psi, phi )
= 0""",
'heat_conduct' :
"""dw_volume_dot.2.Omega( s, dT/dt ) +
dw_laplace.2.Omega( m.thermal_conductivity, s, T )
= dw_electric_source.2.Omega( m.electric_conductivity, s, phi
)""",
'thermo_mech':
"""dw_lin_elastic.2.Omega(m.mech_solid, v, u)
- dw_biot.2.Omega(m.alpha_therm, v, T)
= 0""",
}
#! Solvers
#! -------
solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton', {
'i_max' : 1,
'eps_a' : 1e-10,
'eps_r' : 1.0,
'problem': 'nonlinear',
}),
'ts': ('ts.simple', {
't0': t0,
't1': t1,
'dt': None,
'n_step': n_step,
}),
}
--------------------------------------------------------
Schöne Grüße
kind regards
i.A. Manuel Feuchter
Competence Unit Engineering
-------------------------------------------------------------------------------------------------------------------------------------
Dr.-Ing. Manuel Feuchter
Head of Multibody & Structure Analysis
TWT GmbH - Science & Innovation
Ernsthaldenstraße 17
D-70565 Stuttgart
E-Mail: manuel.feuchter(a)twt-gmbh.de
Mobil: +49 162 / 212 48 16
www.twt-gmbh.de
-------------------------------------------------------------------------------------------------------------------------------------
Geschäftsführung: Dr. Dimitrios Vartziotis, Joachim Laicher (Stv.), Frank
Beutenmüller (Stv.)
Registergericht: Amtsgericht Stuttgart, HRB Nr. 212778
Umsatzsteuer: ID-Nr.: DE147841145
-------------------------------------------------------------------------------------------------------------------------------------
2 years, 1 month
- 2
- 2
Fwd: Re: Re: Several questions for solving a linear elastic problem
by Robert Cimrman
Forwarding to the list...
-------- Forwarded Message --------
Subject: Re: [sfepy] Re: Several questions for solving a linear elastic problem
Date: Tue, 29 Aug 2017 10:59:10 +0200
From: Kathrin Sobe <kathrinsbriefkasten(a)googlemail.com>
To: Robert Cimrman <cimrman3(a)ntc.zcu.cz>
Hi Robert,
just to let you know: This is the result of the example code after I
changed the things according to your suggestions. Thank you.
I think the questions are answered and I also get a nice result picture. In
the next phase I will try the same with the real mesh object.
My modifications are as follows:
* the surface region is extracted directly with:
surface = domain.create_region('Toppart',
'vertices of surface',
'facet') # 528 facets
* the material function returns a force vector for each quadrature point: I
just splitted the triangle force vector by the number of quadrature points
per triangle (facet). I hope this is the right way to distributet the
forces.
* result extraction with post processing: this works after adding
copy_materials = False: The second approach from
examples/linear_elasticity/modal_analysis.py, around line 291 did not work
in my example. The code is there in comment lines. The code was working,
but the result was zero for all elements? However with post_processing I am
able to access the results.
* Furthermore, I changed the definition of the linear elastic material to
solid = Material('solid', values={'D': stiffness_from_youngpoisson(dim=3,
young=200, poisson=0.3)}). In my case the material is defined by the
youngs-modulus and the poissions ratio.
* the gravity is changed to: t3 = Term.new('dw_volume_lvf(f.val, v)',
integral, omega, f=f, v=v) The values for the material are also defined by
a function.
Here is the complete modified code example:
# coding: utf-8
from __future__ import print_function
from __future__ import absolute_import
from argparse import ArgumentParser
import dill
import functools
import numpy
import time
import sys
from sfepy.discrete.common.mappings import get_physical_qps
sys.path.append('.')
from sfepy import data_dir
from sfepy.base.base import IndexedStruct, Struct
from sfepy.discrete import (FieldVariable, Material, Integral, Integrals,
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_youngpoisson
"""u
aus: http://sfepy.org/doc-devel/tutorial.html
"""
help = {
'show' : 'show the results figure',
}
r"""
Linear elasticity with pressure traction load on a surface and constrained
to
one-dimensional motion.
Siehe Beispiel als Skript hier!!!!!!!!!1
http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_trac...
Find :math:`\ul{u}` such that:
.. math::
\int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u})
= - \int_{\Gamma_{right}} \ul{v} \cdot \ull{\sigma} \cdot \ul{n}
\;, \quad \forall \ul{v} \;,
where
.. math::
D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) +
\lambda \ \delta_{ij} \delta_{kl}
\;.
and :math:`\ull{\sigma} \cdot \ul{n} = \bar{p} \ull{I} \cdot \ul{n}`
with given traction pressure :math:`\bar{p}`.
"""
def post_process(out, problem, state, extend=False):
"""
This will be called after the problem is solved.
Parameters
----------
out : dict
The output dictionary, where this function will store additional
data.
problem : Problem instance
The current Problem instance.
state : State instance
The computed state, containing FE coefficients of all the unknown
variables.
extend : bool
The flag indicating whether to extend the output data to the whole
domain. It can be ignored if the problem is solved on the whole
domain already.
Returns
-------
out : dict
The updated output dictionary.
"""
# 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. problem.equations.materials[0]
stress = problem.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)',
mode='el_avg', copy_materials=False)
out['cauchy_stress'] = Struct(name='output_data',
mode='cell', data=stress,
dofs=None) # # size 1782 mal 6 Werte
Pro Element eine Zeile (=1782)
return out
def _get_loadcase():
u"""
Returns a list of facets and a list of force vectors.
The order defines the matching between facet and force vector.
"""
facets = dill.load(open(fname, "r"))
# #force vectors from another example: temporary solution: take the
necessary number of vectors
force_vectors_example = dill.load(open(fname2, "r"))
force_vectors = force_vectors_example[:len(facets)]
return facets, force_vectors
def _set_not_top_region_zero(facets, force_vectors, facet_indices_top):
u"""
Surface facets that are not in the top region get the zero force vector.
"""
zero_vector = numpy.array([0.0, 0.0, 0.0])
for i, f in enumerate(facets):
if f not in facet_indices_top:
force_vectors[i] = zero_vector
else:
pass
return facets, force_vectors
def _get_qp_force_vectors_from_loadcase(force_vectors, facets, region, qps,
number_qps_facet):
u"""
Ensure order of region_facets and divide force vectors for quadrature
points.
Number of quadrature points per facets is always the same? Even if the
triangles have a different size?
"""
force_vectors_qps = numpy.repeat(force_vectors, 3, axis=0)
force_vectors_qps = force_vectors_qps / float(number_qps_facet)
return force_vectors_qps
def linear_tension(ts, coors, mode=None, force_vectors_qp=None, **kwargs):
u"""
mapping between force vectors (one for every triangle) and
quadrature points.
coors are quadrature points: here are three quadrature points per
tetrahedron
that means: the force vector is divided by three and set to every qp
what happens when the number of quadrature points is not even for all
triangles?
"""
if mode == 'qp':
# shape: 128x3
# Array with tension values for each quadrature point (constant
-0.8)
# val = numpy.tile(-10.0, (coor.shape[0], 1, 1)) # numpy.tile(1.0,
(3, 1, 1)): [[[1]], [[1]], [[1]]]
# val = numpy.tile([[0.001], [0.01], [1.8]], (coors.shape[0], 1,
1)) # 1584 mal [[0.5 -0.45284157 0.4803068], ...]
val = numpy.reshape(force_vectors_qp, (coors.shape[0], 3, 1))
return {'val' : val}
def get_gravity(ts, coors, mode=None, **kwargs):
u"""
Gravity in z-direction.
Gravity is a volume load: here constant for every quadrature point.
Mapping also necessary (see above)
"""
if mode == 'qp':
val = numpy.tile([[0.0], [0.0], [10.5]], (coors.shape[0], 1, 1))
return {'val' : val}
def cube_displacement_medium():
u"""
3D Box with tetrahedron elements. Problem type: linear elastic problem
Input:
force vectors for the top triangles (top surface of the tetrahedron
mesh)
and gravity for each tetrahedron.
Output:
- displacement for each vertex
- stress tensor
"""
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()
# load mesh
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)
# why refine?
# domain.refine()
omega = domain.create_region('Omega', 'all') # Gesamtbereich
# facets: bounding box here: -0.5, -0.5, -0.5 bis 0.5, 0.5, 0.5
bb = mesh.get_bounding_box()
top = domain.create_region('Top',
'vertices in (z >= 0.48)',
'facet') # 88 facets
# gives you vertices of each face in the top region, that can be used
to index mesh.coors (mesh.coors[vof]).
vertices_indexes_facets_top = mesh.cmesh.get_incident(0, top.faces,
2).reshape((-1, 3))
# vertices_facets_top = mesh.coors[vertices_indexes_facets_top]
# facet_indices_top = top.get_facet_indices()
surface = domain.create_region('Toppart',
'vertices of surface',
'facet') # 528 facets
# vertices of facets
# defines the order of facets: order has to be applied on facets load
vectors
vertices_indexes_facets_surface = mesh.cmesh.get_incident(0,
surface.faces, 2).reshape((-1, 3))
# gives you vertices of each face in the top region, that can be used
to index mesh.coors (mesh.coors[vof]).
z_min = bb[0][2]
z_max = bb[1][2]
z_range = z_max - z_min
z_tol = 0.01 * z_range
z_min_tol = z_min + z_tol
bottom = domain.create_region('Bottom',
'vertices in (z <= {})'.format(z_min_tol),
'vertex')
field = Field.from_args('displacement', numpy.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)})
# changed to the given matarial parameter e-modulus and poisson
solid = Material('solid', values={'D':
stiffness_from_youngpoisson(dim=3, young=200, poisson=0.3)}) # young=29
oder 200
# # loadcase with forces on the surface
# # currently: 528 elements and 528 force vectors
facets, force_vectors = _get_loadcase()
facets, force_vectors = _set_not_top_region_zero(facets, force_vectors,
vertices_indexes_facets_top)
integral = Integral('i', order=2)
# # quadrature points in the order of surface facets
# # number of quadrature points for each facet defines the splitting
factor?
# # e.g. three qps per facet: divide force vector by three
# # returns a force vector array with one force vector for each qp
qps1 = get_physical_qps(surface, integral)
qps_shape = qps1.shape # (528, 3, 3)
number_qps_facet = qps_shape[1]
force_vectors_qp = _get_qp_force_vectors_from_loadcase(force_vectors,
facets, surface, qps1, number_qps_facet)
force_vectors_qp = force_vectors_qp * numpy.array([0, 0, -1])
linear_tension_function = functools.partial(linear_tension,
force_vectors_qp=force_vectors_qp)
load = Material('load', function=linear_tension_function) # material
defined by values or function
# volume forces f (gravity modelling)
f = Material('f', function=Function('get_gravity', get_gravity))
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, surface,
load=load, v=v)
# old:
# t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid,
v=v, u=u)
t3 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
# ???? is this right??
eq1 = Equation('balance_of_forces', t1 + t2 + t3)
eqs = Equations([eq1])
fixbottom = EssentialBC('fixbottom', bottom, {'u.all' : 0.0})
ls = ScipyDirect({})
nls_status = IndexedStruct()
# configurable residuum
# absolute_tolerance = 1e-10
absolute_tolerance = 0.1
relative_tolerance = 1.0
nls = Newton({ 'i_max' : 1,
'eps_a' : absolute_tolerance,
'eps_r' : relative_tolerance,
'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.set_materials([load, solid])
pb.save_regions_as_groups('regions')
# ebcs: elastic boundary conditions
pb.time_update(ebcs=Conditions([fixbottom]))
start_time = time.time()
vec = pb.solve()
run_time = time.time()-start_time
print("runtime medium problem {} {}: {}".format(mesh.n_nod, mesh.n_el,
run_time))
print(nls_status)
# von mises stress berechnen: is working now!!! copy_materials = False
in function!
# pb.save_state('linear_elastic_traction_cylinder.vtk', vec,
post_process_hook=post_process)
name = 'linear_elastic_traction_medium_box.vtk'
pb.save_state(name, vec, post_process_hook=post_process)
######### stress and strain results are zero #######
# # alternative for calculating the stress tensor
# state = pb.create_state()
# aux = state.create_output_dict()
# strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)',
# mode='el_avg', verbose=False)
#
#
# aux['cauchy_strain'] = Struct(name='output_data',
# mode='cell', data=strain,
# dofs=None)
#
#
# stress = pb.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)',
# mode='el_avg', verbose=False,
copy_materials=False)
#
# aux['cauchy_stress'] = Struct(name='output_data',
# mode='cell', data=stress,
# dofs=None) # # size 1782 mal 6 Werte
Pro Element eine Zeile (=1782)
if options.show:
view = Viewer(name)
view(vector_mode='warp_norm', rel_scaling=2,
is_scalar_bar=True, is_wireframe=True)
pb.reset()
u.reset()
v.reset()
return run_time, mesh.n_nod
if __name__ == '__main__':
runtime_small_cube, nodes_small_cube = cube_displacement_medium()
2017-08-25 17:03 GMT+02:00 Robert Cimrman <cimrman3(a)ntc.zcu.cz>:
> Hi Kathrin,
>
>
> On 08/25/2017 04:42 PM, Kathrin Sobe wrote:
>
>> Hi Robert,
>>
>> thanks a lot for your help.
>>
>>
>> 2017-08-25 14:37 GMT+02:00 Robert Cimrman <cimrman3(a)ntc.zcu.cz>:
>>
>> Hi Kathrin,
>>>
>>> On 08/25/2017 10:37 AM, Kathrin Sobe wrote:
>>>
>>> Hi Robert,
>>>>
>>>> I try to solve a linear elastic problem. I give a example in a python
>>>> file that is a modified version of the example given in
>>>> http://sfepy.org/doc-devel/examples/linear_elasticity/linear
>>>> _elastic_tractions.html
>>>>
>>>> The problem I want to solve is the following:
>>>> I have a object, defined as a 3D mesh of tetrahedrons. In the test
>>>> example I take the medium size 3D cube. The real 3D object is much
>>>> bigger
>>>> and has a more structured surface.
>>>>
>>>> On the top of the object there are certain forces. These forces are
>>>> defined as 3D vectors, one for each surface triangle. Additionally the
>>>> object has certain material parameters given by density, poisson ratio
>>>> and
>>>> elastic modulus.
>>>>
>>>> Working through the linear elastic examples I figured out two different
>>>> types of solutions:
>>>> 1. modeling the forces at the surface as a material with a function
>>>> (manual 4.3.9)
>>>> 2. modeling the forces at the surface as a boundary condition
>>>>
>>>> So the first question is in general: Is this the right way to model my
>>>> problem? At the moment I think solution type number one is the better
>>>> approach?
>>>>
>>>>
>>> Yes. The dw_surface_ltr term you use corresponds to the Neumann boundary
>>> conditions = applied forces/tractions/pressure in the weak form of the
>>> elasticity PDE. Because you need a position-dependent value of the load,
>>> you need to use a material parameter given by a function.
>>>
>>
>>
>> ok, than I will concentrate on that approach. Additionally, I want to
>> model
>> the gravity. For that I added the rho parameter to the material defintion
>> and defined a mass equation like this:
>> # # rho ist die Dichte: bei geringerer Dichte, wird der Block mehr
>> verformt: ist klar
>> solid = Material('solid', values={'D': stiffness_from_lame(dim=3,
>> lam=5.769, mu=3.846), 'rho': 0.501})
>> load = Material('load', function=linear_tension) # material defined
>> by
>> values or function
>>
>> 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, surface,
>> load=load, v=v)
>> t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid,
>> v=v, u=u)
>> eq1 = Equation('elasticity', t1 - t2)
>> eq2 = Equation('mass', t3)
>> eqs = Equations([eq1, eq2])
>>
>> I am not sure if this is the right way to add the gravity to the problem.
>>
>
> Use the dw_volume_lvf term for volume forces just as you use
> dw_surface_ltr for surface forces. The material parameter here is a force
> per unit volume. The mass matrix would be useful in a dynamic simulation,
> not here.
>
>
>>> Going through the first approach I have the following more detailed
>>>
>>>> problems:
>>>> * Region definition: I need a way to define the surface triangles of the
>>>> 3D object. Because of the structure of the object, this cannot be done
>>>> with
>>>> coordinate ranges.
>>>>
>>>>
>>> Is it going to be applied on the whole surface? Anyway, you can select
>>> the
>>> whole surface of a mesh with 'vertices of surface', for example:
>>>
>>> surf = domain.create_region('surf', 'vertices of surface', 'facet')
>>>
>>> Check the manual section 4.3.3 (on regions) for additional possibilities
>>> of defining regions and combining existing regions using various
>>> operations.
>>>
>>> Or, if your mesh comes from some mesh generator/cad software, you can
>>> mark
>>> various vertex / cell groups there and use those for region definitions.
>>> This depends on the mesher/format you use - let me know if you need help
>>> with that.
>>>
>>>
>> The surface extraction is exactly what I need, cool. I found that out as
>> well. Thanks for clarifying the indexing. I will try that.
>>
>> Just as an information: Not all surface triangles have a force vector,
>> only
>> the triangles in the upper region of the geometry.
>>
>
> Then you may need to remove the extra faces from the surface region (or
> simply define the force zero there).
>
>
>>> * Material definition with force vectors: the 3D force vectors are given
>>>
>>>> for each surface triangle (defined at the centroid of each triangle). In
>>>> the material function examples I found only the definition of forces at
>>>> quadrature points and only as scalar values.
>>>>
>>>>
>>> You can use vectors there too, for example:
>>>
>>> val = nm.tile([[-0.8], [0.1], [0.2]], (coor.shape[0], 1, 1))
>>>
>>> To have a constant value per element, you need to know in which element
>>> or
>>> face each quadrature point is. This information is available in the
>>> material function (for example, kwargs['term'].region.faces contains the
>>> list of faces of the region, etc.), but the actual code depends on how
>>> your
>>> material function is represented - is it a table of (cell-value) pairs,
>>> or
>>> something else?
>>>
>>>
>> At the moment my force defintion is exacly as you expected: cell-value
>> pairs. The cell is the surface triangle and the value is the force vector.
>>
>> Is it that easy: Find the cell for each quadrature point and return the
>> cell force vector as its load? Do I need to split the force vector of the
>> cell between its quadrature points? Is there always exacly one cell for
>> each quadrature point?
>>
>
> You need to provide a force per unit area vector (so yes, we may call it
> splitting) in each quadrature point. Each quadrature point is in exactly
> one face, yes. The number of quadrature point per face can be found as
> follows:
>
> qps = t2.get_physical_qps()
> print qps.shape
> (88, 3, 3)
>
> = 88 faces, 3 QPs per face, 3 coordinates (in 3D)
>
> So this tells you, that the coor argument of linear_tension has 3 rows for
> each face (in the order of top.faces (= t2.region.faces).
>
> r.
>
> _______________________________________________
> SfePy mailing list
> sfepy(a)python.org
> https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/
>
2 years, 1 month
- 1
- 1
Re: Several questions for solving a linear elastic problem
by Robert Cimrman
Hi Kathrin,
On 08/25/2017 04:42 PM, Kathrin Sobe wrote:
> Hi Robert,
>
> thanks a lot for your help.
>
>
> 2017-08-25 14:37 GMT+02:00 Robert Cimrman <cimrman3(a)ntc.zcu.cz>:
>
>> Hi Kathrin,
>>
>> On 08/25/2017 10:37 AM, Kathrin Sobe wrote:
>>
>>> Hi Robert,
>>>
>>> I try to solve a linear elastic problem. I give a example in a python
>>> file that is a modified version of the example given in
>>> http://sfepy.org/doc-devel/examples/linear_elasticity/linear
>>> _elastic_tractions.html
>>>
>>> The problem I want to solve is the following:
>>> I have a object, defined as a 3D mesh of tetrahedrons. In the test
>>> example I take the medium size 3D cube. The real 3D object is much bigger
>>> and has a more structured surface.
>>>
>>> On the top of the object there are certain forces. These forces are
>>> defined as 3D vectors, one for each surface triangle. Additionally the
>>> object has certain material parameters given by density, poisson ratio and
>>> elastic modulus.
>>>
>>> Working through the linear elastic examples I figured out two different
>>> types of solutions:
>>> 1. modeling the forces at the surface as a material with a function
>>> (manual 4.3.9)
>>> 2. modeling the forces at the surface as a boundary condition
>>>
>>> So the first question is in general: Is this the right way to model my
>>> problem? At the moment I think solution type number one is the better
>>> approach?
>>>
>>
>> Yes. The dw_surface_ltr term you use corresponds to the Neumann boundary
>> conditions = applied forces/tractions/pressure in the weak form of the
>> elasticity PDE. Because you need a position-dependent value of the load,
>> you need to use a material parameter given by a function.
>
>
> ok, than I will concentrate on that approach. Additionally, I want to model
> the gravity. For that I added the rho parameter to the material defintion
> and defined a mass equation like this:
> # # rho ist die Dichte: bei geringerer Dichte, wird der Block mehr
> verformt: ist klar
> solid = Material('solid', values={'D': stiffness_from_lame(dim=3,
> lam=5.769, mu=3.846), 'rho': 0.501})
> load = Material('load', function=linear_tension) # material defined by
> values or function
>
> 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, surface,
> load=load, v=v)
> t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid,
> v=v, u=u)
> eq1 = Equation('elasticity', t1 - t2)
> eq2 = Equation('mass', t3)
> eqs = Equations([eq1, eq2])
>
> I am not sure if this is the right way to add the gravity to the problem.
Use the dw_volume_lvf term for volume forces just as you use dw_surface_ltr for
surface forces. The material parameter here is a force per unit volume. The
mass matrix would be useful in a dynamic simulation, not here.
>>
>> Going through the first approach I have the following more detailed
>>> problems:
>>> * Region definition: I need a way to define the surface triangles of the
>>> 3D object. Because of the structure of the object, this cannot be done with
>>> coordinate ranges.
>>>
>>
>> Is it going to be applied on the whole surface? Anyway, you can select the
>> whole surface of a mesh with 'vertices of surface', for example:
>>
>> surf = domain.create_region('surf', 'vertices of surface', 'facet')
>>
>> Check the manual section 4.3.3 (on regions) for additional possibilities
>> of defining regions and combining existing regions using various operations.
>>
>> Or, if your mesh comes from some mesh generator/cad software, you can mark
>> various vertex / cell groups there and use those for region definitions.
>> This depends on the mesher/format you use - let me know if you need help
>> with that.
>>
>
> The surface extraction is exactly what I need, cool. I found that out as
> well. Thanks for clarifying the indexing. I will try that.
>
> Just as an information: Not all surface triangles have a force vector, only
> the triangles in the upper region of the geometry.
Then you may need to remove the extra faces from the surface region (or simply
define the force zero there).
>>
>> * Material definition with force vectors: the 3D force vectors are given
>>> for each surface triangle (defined at the centroid of each triangle). In
>>> the material function examples I found only the definition of forces at
>>> quadrature points and only as scalar values.
>>>
>>
>> You can use vectors there too, for example:
>>
>> val = nm.tile([[-0.8], [0.1], [0.2]], (coor.shape[0], 1, 1))
>>
>> To have a constant value per element, you need to know in which element or
>> face each quadrature point is. This information is available in the
>> material function (for example, kwargs['term'].region.faces contains the
>> list of faces of the region, etc.), but the actual code depends on how your
>> material function is represented - is it a table of (cell-value) pairs, or
>> something else?
>>
>
> At the moment my force defintion is exacly as you expected: cell-value
> pairs. The cell is the surface triangle and the value is the force vector.
>
> Is it that easy: Find the cell for each quadrature point and return the
> cell force vector as its load? Do I need to split the force vector of the
> cell between its quadrature points? Is there always exacly one cell for
> each quadrature point?
You need to provide a force per unit area vector (so yes, we may call it
splitting) in each quadrature point. Each quadrature point is in exactly one
face, yes. The number of quadrature point per face can be found as follows:
qps = t2.get_physical_qps()
print qps.shape
(88, 3, 3)
= 88 faces, 3 QPs per face, 3 coordinates (in 3D)
So this tells you, that the coor argument of linear_tension has 3 rows for each
face (in the order of top.faces (= t2.region.faces).
r.
2 years, 1 month
- 1
- 0
Several questions for solving a linear elastic problem
by Kathrin Sobe
Hi Robert,
I try to solve a linear elastic problem. I give a example in a python file that is a modified version of the example given in http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_trac...
The problem I want to solve is the following:
I have a object, defined as a 3D mesh of tetrahedrons. In the test example I take the medium size 3D cube. The real 3D object is much bigger and has a more structured surface.
On the top of the object there are certain forces. These forces are defined as 3D vectors, one for each surface triangle. Additionally the object has certain material parameters given by density, poisson ratio and elastic modulus.
Working through the linear elastic examples I figured out two different types of solutions:
1. modeling the forces at the surface as a material with a function (manual 4.3.9)
2. modeling the forces at the surface as a boundary condition
So the first question is in general: Is this the right way to model my problem? At the moment I think solution type number one is the better approach?
Going through the first approach I have the following more detailed problems:
* Region definition: I need a way to define the surface triangles of the 3D object. Because of the structure of the object, this cannot be done with coordinate ranges.
* Material definition with force vectors: the 3D force vectors are given for each surface triangle (defined at the centroid of each triangle). In the material function examples I found only the definition of forces at quadrature points and only as scalar values.
* After solving the problem I need to extract the node displacement and the stress. For the stress calculation I tried to use a post_process function as visible in the example. But this was not working.
Thank you for your help. I hope the code in my example code makes the problems more clear.
Regards,
Kathrin
******************** code example: solution type with force vectors defined in material parameter *******************
# 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
"""u
aus: http://sfepy.org/doc-devel/tutorial.html
"""
help = {
'show' : 'show the results figure',
}
r"""
Linear elasticity with pressure traction load on a surface and constrained to
one-dimensional motion.
Siehe Beispiel als Skript hier!!!!!!!!!1
http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_trac...
Find :math:`\ul{u}` such that:
.. math::
\int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u})
= - \int_{\Gamma_{right}} \ul{v} \cdot \ull{\sigma} \cdot \ul{n}
\;, \quad \forall \ul{v} \;,
where
.. math::
D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) +
\lambda \ \delta_{ij} \delta_{kl}
\;.
and :math:`\ull{\sigma} \cdot \ul{n} = \bar{p} \ull{I} \cdot \ul{n}`
with given traction pressure :math:`\bar{p}`.
"""
def post_process(out, problem, state, extend=False):
"""
This will be called after the problem is solved.
Parameters
----------
out : dict
The output dictionary, where this function will store additional
data.
problem : Problem instance
The current Problem instance.
state : State instance
The computed state, containing FE coefficients of all the unknown
variables.
extend : bool
The flag indicating whether to extend the output data to the whole
domain. It can be ignored if the problem is solved on the whole
domain already.
Returns
-------
out : dict
The updated output dictionary.
"""
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
# # coor: 264 x3 --> qp points in top region!
def linear_tension(ts, coor, mode=None, **kwargs):
if mode == 'qp':
# shape: 128x3
# Array with tension values for each quadrature point (constant -0.8)
val = nm.tile(-0.8, (coor.shape[0], 1, 1)) # numpy.tile(1.0, (3, 1, 1)): [[[1]], [[1]], [[1]]]
# ## with vectors it does not work
# val = nm.tile(nm.array([-0.01, 0.05, -0.8]), (coor.shape[0], 1, 1))
return {'val' : val}
def cube_displacement_medium():
u"""
3D Box with tetrahedron elements. Problem type: linear elastic problem
linear element approach (later quadratic)
Input:
I have force vectors for the top triangles (top surface of the tetrahedron mesh)
and gravity for each tetrahedron.
Output:
- displacement for each vertex
- stress tensor
"""
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()
# load mesh
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 mit Mesh definieren: hier 448 Knoten, 1782 Elemente
domain = FEDomain('domain', mesh)
# why refine?
domain.refine()
omega = domain.create_region('Omega', 'all') # Gesamtbereich
# facets: bounding box here: -0.5, -0.5, -0.5 bis 0.5, 0.5, 0.5
bb = mesh.get_bounding_box()
top = domain.create_region('Top',
'vertices in (z >= 0.48)',
'facet')
bottom = domain.create_region('Bottom',
'vertices in (z <= -0.48)',
'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) # material defined by values or function
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, top, load=load, v=v)
eq = Equation('elasticity', t1 - t2)
eqs = Equations([eq])
fixbottom = EssentialBC('fixbottom', bottom, {'u.all' : 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')
# ebcs: elastic boundary conditions
pb.time_update(ebcs=Conditions([fixbottom]))
start_time = time.time()
vec = pb.solve()
run_time = time.time()-start_time
print("runtime medium problem {} {}: {}".format(mesh.n_nod, mesh.n_el, run_time))
print(nls_status)
# von mises stress calculation (not working)
pb.save_state('linear_elastic_simulation2.vtk', vec)
# pb.save_state('linear_elastic_traction_cylinder.vtk', vec, post_process_hook=post_process)
if options.show:
view = Viewer('linear_elastic_simulation2.vtk')
view(vector_mode='warp_norm', rel_scaling=2,
is_scalar_bar=True, is_wireframe=True)
pb.reset()
u.reset()
v.reset()
return run_time, mesh.n_nod
if __name__ == '__main__':
runtime_small_cube, nodes_small_cube = cube_displacement_medium()
2 years, 1 month
- 2
- 3