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),
}
Dear,
Recently I have used SfePy for the computational homogenization and I found
that SfePy has offered the homogenization component. Thanks for that.
However, this feature is only for microscopic structure that means the
simulation is for an RVE (Representative Volume Elements) with the
deformation gradient or strain in the context of solid mechanics as the
input to set boundary conditions for microscopic structure.
In my work, I simulate a structure made of the inhomogeneous materials.
This means that at each Gauss points in computing of stiffness matrix we
compute the average stress and effective moduli by simulating the
microscopic RVE (this feature is already available in SfePy). However, In
SfePy we set up a material object with its parameters at the beginning and
call or setup a *terms *file.
Do you know the solution or the way using terms feature to compute the
homogenization in macro scale?
Thanks and Best Regards,
Minh Nguyen
Hi,
When I ran the time_poisson_interactive [1] example, I noticed that probing took a lot of time. Speciffically, when I:
- used pb.get_tss_functions(state0, save_results=False) to turn off saving to disk,
- commented out any matplotlib plotting part, and
- turned on options.probe
- ran 201 time steps
... the problem was solved in 46.9 seconds (the tss() call at the bottom), but probing itself (only the entire for loop for the two probes from "for ii, probe in enumerate(probes):" to "append(results)") took 27.1 seconds in total (0.135 sec each step on average).
Further profiling reveals that everytime when Field.evalute_at() is called during probing, the get_ref_coors() function and its subsequent downstream functions are called to (what I understand) find the evalutaion points in the mesh, taking up 26.4 seconds in total.
This is rather inefficient. Is there a way to reuse the previously calculated probing cells, since the probes and mesh are constant throughout the problem?
[1] http://sfepy.org/doc-devel/examples/diffusion/time_poisson_interactive.html
Regards,
Bill
Hi,
I tried the compare_elastic_materials.py as below and got errors.
$ ./simple.py examples/large_deformation/compare_elastic_materials.py
sfepy: left over: ['mesh_hook', 'stiffness_from_lame', 'UserMeshIO', 'gen_block_mesh', 'verbose', '_filename']
sfepy: reading mesh [user] (function:mesh_hook)...
sfepy: ...done in 0.00 s
sfepy: creating regions...
sfepy: Omega
sfepy: Bottom
sfepy: Top
sfepy: ...done in 0.04 s
sfepy: equation "linear":
sfepy: dw_lin_elastic.i.Omega(solid.D, v, u)
= dw_surface_ltr.isurf.Top(load.val, v)
sfepy: equation "neo-Hookean":
sfepy: dw_tl_he_neohook.i.Omega(solid.mu_nh, v, u)
+ dw_tl_bulk_penalty.i.Omega(solid.K, v, u)
= dw_surface_ltr.isurf.Top(load.val, v)
sfepy: equation "Mooney-Rivlin":
sfepy: dw_tl_he_neohook.i.Omega(solid.mu_mr, v, u)
+ dw_tl_he_mooney_rivlin.i.Omega(solid.kappa, v, u)
+ dw_tl_bulk_penalty.i.Omega(solid.K, v, u)
= dw_surface_ltr.isurf.Top(load.val, v)
sfepy: using solvers:
ts: ts
nls: newton
ls: ls
sfepy: updating variables...
sfepy: ...done
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: matrix shape: (28, 28)
sfepy: assembling matrix graph...
sfepy: ...done in 0.00 s
sfepy: matrix structural nonzeros: 688 (8.78e-01% fill)
sfepy: updating materials...
sfepy: load
Traceback (most recent call last):
File "./simple.py", line 177, in <module>
main()
File "./simple.py", line 174, in main
app()
File "/home/pyontaku14/anaconda3/sfepy/sfepy/applications/application.py", line 29, in call_basic
return self.call(**kwargs)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/applications/pde_solver_app.py", line 218, in call
post_process_hook_final=self.post_process_hook_final)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/problem.py", line 1383, in solve
status=status)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/solvers/ts_solvers.py", line 34, in _standard_ts_call
status=status, **kwargs)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/solvers/ts_solvers.py", line 142, in __call__
vec0 = init_fun(ts, vec0)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/problem.py", line 1187, in init_fun
self.init_time(ts)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/problem.py", line 765, in init_time
self.update_materials(mode='force')
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/problem.py", line 575, in update_materials
problem=self, verbose=verbose)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/equations.py", line 340, in time_update_materials
verbose=verbose)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/materials.py", line 59, in time_update
mat.time_update(ts, equations, mode=mode, problem=problem)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/materials.py", line 304, in time_update
self.update_data(key, ts, equations, term, problem=problem)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/materials.py", line 221, in update_data
**self.extra_args)
File "/home/pyontaku14/anaconda3/sfepy/sfepy/discrete/functions.py", line 36, in __call__
return self.function(*args, **_kwargs)
TypeError: <lambda>() got an unexpected keyword argument 'equations'
$
The version is the newest.
$ ./simple.py --version
simple.py 2018.1+git.4da6f2a4a87c44869ad7016686ccddc1e75bb845
Could you please give me advice to solve the error?
Best regards,
Takuo Fujita
Hi Robert,
I saw this conversation with Santiago and I am having what I think a very similar problem.
I would like to define Material property based on the defined regions in meshes. In the examples I saw, the material property is always defined either as a constant or based on the coordinates stored in "coors".
For example, in the poisson with source term example[1], the load is defined as "5e5 * y coordinate". How can I define the load as say 5.5e5 for vertices inside the 'Omega_L' region, and zero outside of Omega_L, thus effectively creating a constant heat source in the middle ball? In other words, how can I get the point indices in "coors" that belong to the "Omega_L" region?
I understand that "Omega_L" in this example is defined by the "get_middle_ball" function so one can just implement that function inside the material property defining function to find which "coors" are inside the ball. However for my purpose my sources are defined as groups in the mesh, where I can retrieve the region by doing "domain.create_region(name='source_region', select='cells of group 301')".
This would be an example mesh that I would generate, where cells of group 301 correspond to a heat source with a finite heat source Q, and the rest of the space does not generate heat:
MeshVersionFormatted 2
Dimension
3
Vertices
28
0 0 0.5 1
0 0 0 2
0 1 0.5 3
###More Vertices###
Tetrahedra
48
4 8 20 21 301
18 8 6 21 301
###More tetrahedra labelled as group 301, this is a heat source with material property A###
5 18 6 19 301
11 9 23 28 302
14 16 27 24 302
###More tetrahedra, these have material property B###
25 13 14 24 302
End
[1] http://sfepy.org/doc-devel/examples/diffusion/poisson_functions.html
Thank you for your help!
Regards,
Bill
Hi,
I am a new user of SfePy and finite element programming in general, and I
have started my own finite element project last week. After comparing SfePy
with other packages like FEniCS and GetFEM++ I have decided to give SfePy a
try. My problem is close to a diffusion problem. I have a fair amount of
scientific Python experience in general, including some VTK.
I know I can export the computation result to VTK files in Unstructured
Grid format, which is performed at the end of many examples. My question
is: how can I access the final result, at each vertex, directly from the
Python program itself, saving the result into say a Numpy array or
equivalent data structure.
For example, the Time Thermal Interactive example
<http://sfepy.org/doc-devel/examples/diffusion/time_poisson_interactive.html>
exports domain.01.vtk to domain.02.vtk, containing points, cells, and
temperature at each point ("SCALARS T float 1.") Is there a way to access
these numbers directly at the end of the computation, from say the Problem
variable pb or State variable state0? I feel like this is probably an easy
task but I could not figure out on my own just by looking into the State
and Problem variables.
I need these values directly in Python as I have other part of my model
that takes the temperature values at these vertices. Since the other
computation will use the same mesh I would like to have the mesh and all
solved values.
Thanks!
Bill
