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, 1 month

1D line elements in SfePy
by Nimish

I am currrently looking for FEM packages to help me solve a system of
beams and columns, basically a collection of 1D bernoulli/timoshenko
line elements.
I started reading SfePy docs and i am getting the idea that doing the
above is not really possible here, am i right?
Are only 2D area elements permitted in SfePy?
Or is there any direct support for solving 1D line elements too..
Cheers
Nimish
3 years, 2 months

Evaluate a solution in the arbitrary point in the domain
by Alec Kalinin

Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node, but in
any arbitrary point in the domain with the given (x, y, z) coordinates?
For example, consider Dirichlet problem for Poisson equation. We apply
essential boundary conditions on the surface nodes and after the problem
has been solved we have the solution vector, i.e. vector of values in the
FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is
not FEM mesh node. What is the best way to obtain solution in this point v?
Sincerely,
Alec Kalinin
5 years

elliptical hole under isotropic stress
by David N. Mashburn

Hello sfepy developers and users!
I am modelling a simple linear elastic sheet under isotropic stress with
an elliptical hole in the center (and I have it working under sfepy,
great little platform!).
It is obvious the model should initially yield more easily in the
direction of the short axis of the ellipse. What is not so obvious to me
is what should happen in the limit as stress goes to infinity. Part of
me wants to believe that the hole should eventually become a circular,
but the results of the simulation show that the ellipse eventually
switches its aspect ratio with what was the the short axis becoming the
long axis and vice-versa.
My question is whether:
A: The finite element result is the product of a
small-displacement/non-moving mesh artifact (and if so, if there is a
way to get the correct behavior using sfepy...)
OR
B: My intuition about the physical behavior of this ideal system is
incorrect and the ellipse really wouldn't round out into a circle under
increasingly large stress (aka, the FE model is still physical/correct
with large displacements).
This might be obvious to people who have done more finite element
modeling than I have, but thanks anyway! I'm attaching a picture to make
it easier to see at a glance (quarter-ellipse with x and y symmetry
boundary conditions and equal tractions applied at the top and right
boundaries).
Thanks!
-David Mashburn
5 years, 3 months

ANN: SfePy 2013.4
by Robert Cimrman

I am pleased to announce release 2013.4 of SfePy.
Description
-----------
SfePy (simple finite elements in Python) is a software for solving
systems of coupled partial differential equations by the finite element
method. The code is based on NumPy and SciPy packages. It is distributed
under the new BSD license.
Home page: http://sfepy.org
Mailing list: http://groups.google.com/group/sfepy-devel
Git (source) repository, issue tracker, wiki: http://github.com/sfepy
Highlights of this release
--------------------------
- simplified quadrature definition
- equation sequence solver
- initial support for 'plate' integration/connectivity type
- script for visualization of quadrature points and weights
For full release notes see http://docs.sfepy.org/doc/release_notes.html#id1
(rather long and technical).
Best regards,
Robert Cimrman and Contributors (*)
(*) Contributors to this release (alphabetical order):
Vladimír Lukeš, Jaroslav Vondřejc
5 years, 4 months

Re: [sfepy-devel] Re: post process energy calculation
by Robert Cimrman

On 11/15/2013 10:13 AM, Anastasia Konovalova wrote:
> Thanks for the answer, Robert! And excuse me for the long response.
>
> I'm not sure, that I understand you clearly.
> I have added changes from Git, and the computation of determinant in my
> programm is correct now. And your "black magic" code don't produce error
> any more. So you repaired this old bag, right?) And strain can be used for
> further computations now? Thank you very much!
Yes, it should work ok now.
> As for the terms and issues, I compute determinant and other values for the
> energy calculation. If it is possible to make terms for the energy density
> function for hyperelastic material, it would be great and I'll add the
> issue as soon as I can:)
Yes, that would be possible. Add the issue, please, so that it is not forgotten.
Cheers,
r.
> четверг, 7 ноября 2013 г., 16:34:59 UTC+6 пользователь Robert Cimrman
> написал:
>>
>> Hi Anastasia,
>>
>> great job, you have just found a very old bug!
>>
>> The following code (to be added in stress_strain() after i3 computation)
>> should
>> not raise errors but it does:
>>
>> # black magic here...
>> cache =
>> state.variables['u'].evaluate_cache['tl_common'][0][('__v_o1', 0,
>> 0, 'volume', None)]
>>
>> from sfepy.linalg import dot_sequences as dot
>> from sfepy.mechanics.tensors import get_full_indices
>>
>> c = dot(cache.mtx_f, cache.mtx_f, 'ATB')
>> e = 0.5 * (c - nm.eye(3))
>>
>> ii = get_full_indices(3)
>> c2 = cache.sym_c[..., ii, 0]
>> e2 = cache.green_strain[..., ii, 0]
>>
>> assert(nm.allclose(nm.sqrt(i3), cache.det_f, rtol=1e-15, atol=1e-15))
>> assert(nm.allclose(c, c2, rtol=1e-15, atol=1e-15))
>> assert(nm.allclose(e, e2, rtol=1e-15, atol=1e-15))
>>
>> The last line assertion fails, so the Green strain e2 as computed by the
>> dw_tl_he_neohook term call differs from e computed by definition (0.5 *
>> (F^T F
>> - I)).
>>
>> It has escaped me for so long as we have used the green strain only in
>> postprocessing to make "colorful images"... The Green strain has not been
>> used
>> in any further computations. The stresses are ok - they are computed from
>> the C
>> invariants directly.
>>
>> The problem (doubled out-of diagonal entries of E) is now fixed in the
>> repo and
>> your example gives positive i3 in all time steps.
>>
>> Thanks again for noticing the problem!
>>
>> Cheers,
>> r.
>>
>> PS: If you would like to have a term/terms for all the attributes of the
>> evaluate cache above (do 'print cache' to see them), add an issue, please.
>>
>> On 11/06/2013 08:45 AM, Anastasia Konovalova wrote:
>>> Hello Robert!
>>> I got a question about the energy calculation again :)
>>>
>>> If you remember, in my programm there is function get_Wh that is called
>> in
>>> postproc. This function compute the determinant of the right
>> Cauchy-Green
>>> deformation tensor. When I was trying to do a large deformation, a
>>> situation arose, where in some cells the determinant is very small or
>> even
>>> negative.
>>>
>>> I can analyze the results on these elements in Paraview, so I see that
>> the
>>> field of displacements of these cells are normal. I can deform my sample
>> by
>>> displacement, so I find out that the volume of element is normal too,
>>> determinant haven't been far from 1. But If I calculate the determinant
>>> from the components of strain tensor by hands, I'll get this negative
>>> value, so that the computation of the determinant is right too.
>>>
>>> It seems to me that this is the problem of the Green strain tensor
>>> computation.
>>>
>>> I attached the screenshot from Paraview, where you can see the
>> determinant
>>> value, the volume of element, and components of the Green strain tensor,
>>> the file of programm and a simple mesh, and file with the results (for
>>> paraview or some other similar programm) of one of the first steps, when
>>> the negative determinant appears.
>>>
>>> At the function get_Wh I get the the Green strain tensor, compute the
>> right
>>> Cauchy-Green deformation tensor (C = 2E+I), and get i3 =
>> nm.linalg.det(C).
>>> I need the sqrt of i3 to get a J = dV/dV_0, and the error appears there.
>>>
>>> Perhaps, it's too difficult to understand all this, so I will be glad to
>>> any advices. Thank you anyway.
>>>
>>
>>
>
5 years, 5 months

post process energy calculation
by Anastasia Konovalova

Hello, sfepy group,
I have a question about post-process calculations. For my problem I used
the modified examples linear_elastic.py and hyperelastic.py. After the
calculation of strain and stress tensors I need to get energy of
deformation
\int_{V} S : C,
where C = 2*E+I.
I think I need to calculate the multiplication of tensors first, as some
material variable, and then integrate it.
So I have 2 questions:
1. How can I calculate the double scalar multiplication of two arrays
strain and stress that have shape (number of elements, 1, 6, 1)? (Double
scalar multiplication of A and B is a scalar that can be got as
A_{ij}B_{ji}). Because of the shape of arrays I can't use
numpy.inner(stress,strain). I made some function for these calculation, but
perhaps, the simplier way exsists.
2. Then I get the multiplication - the scalar function of energy I should
integrate it. For this I get the dw_volume_integrate term. I'm trying to
integrate material parameter, where energy has been saved.
solid.update({
'E' : get_W(stress,strain),#energy in elements
})
U = problem.evaluate('dw_volume_integrate.i1.Omega(solid.E, v)',)
Or may be I should use some new variable and initialize it in post-proc.
I'm asking for help in correct formulation of this integral. Should I
define energy as material parameter or some new variable? And how I can do
it.
Thanks in advance,
5 years, 5 months

simplified integral definition
by Robert Cimrman

Hi,
the code has just been updated to allow the following definitions of integrals:
integral_1 = {
'name' : 'i1',
'order' : 2,
}
integrals = {
'i1' : 2,
}
There are no more 'kind' keyword, and 'kind' attribute of Integral. They were
removed because the integration kind is fully determined by the region topology.
The old versions
integral_1 = {
'name' : 'i1',
'kind' : 'v',
'order' : 2,
}
integrals = {
'i1' : ('v', 2),
}
still work for backward compatibility, but the kind is ignored and I recommend
to remove it from your files.
Cheers,
r.
5 years, 5 months