On Thursday, July 31, 2014 8:49:23 PM UTC+8, Robert Cimrman wrote:
On 07/31/2014 02:27 PM, Ouyang wrote:
On Thursday, July 31, 2014 6:56:46 PM UTC+8, Robert Cimrman wrote:
Hi Robert,
1. how to output principal stresses in sfepy? I can calculate
On 07/31/2014 12:12 PM, Ouyang wrote: principal
stresses based on \sigma_{ij} in svar below. But I guess it is not straightforward for the task in my question 2. ... sfield = Field.from_args('stress', nm.float64, (6,), pb.domain.regions['Omega']) svar = FieldVariable('sigma', 'parameter', sfield,primary_var_name='(set-to-None)') svar.set_data_from_qp(stress, ivn) ...
I guess you do not need to interpolate the stresses into a field
variable
- you can work with the element-averaged values, no?
I want to evaluate the stresses at center of the block which point load is applied. For brittle materials, rock, we concern its tensile stress (normally minimum stress defined in rock mechanics) at the block center to see if it reaches its tensile strength. So, principal stresses is normally easy to be used in a failure criteria in engineering.
I see. So set_data_from_qp() might be good for you. What it actually does is that it integrates a quantity (given in qp) in all elements around a vertex and divides it by the volume of those elements. So it's a volume-weighted average.
I guess the current element types available in sfepy are constant-strain or constant-stress in a element. So my understanding is that the constant-strain or constant-stress is element-averaged values.
Not only constant-strain - this depends on the field approximation order. If you set the displacement field order to something > 1, even for tetrahedra the train in an element will not be constant.
BTW. a function for computing principal stresses would be a nice
addition
to sfepy.mechanics.tensors.
I am not sure I am good enough for this addition. I can try it when I have time.
Thanks! (Are you interested just in principal stress values, or also in the directions?)
Yes. both magnitudes and directions. I used to use "eigenvalues and eigenvectors" http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#eigenvalues-a... from scipy:
import numpy as np from scipy import linalg A = np.array([[1,5,2],[2,4,1],[3,6,2]]) la,v = linalg.eig(A)
la are the magnitudes of principal stresses and v holds its directions.
Does this "fit" somewhere in tensors.py of sfepy?
I guess not that easy.
2. are there any function to output just a single component of u_i
(displacements) or \sigma_{ij} into vtk file?
Not per se, but it's really easy - just add that to the output dictionary:
from sfepy.mechanics.tensors import get_full_indices s2f = nm.array(get_full_indices(3)) ii = s2f[0, 1] out['cauchy_stress_01'] = Struct(name='output_data', mode='cell', data=stress[..., ii:ii+1, :]) out['u_0'] = Struct(name='output_data', mode='vertex', data = out['u'].data[:, 0:1])
Thanks for this. I will test it.
3. how to output global coordinates for \sigma_{ij} in svar at qp?
You mean how to get the physical coordinates of all the qp in the whole domain? Or something else?
Yes, the physical coordinates associated with the stresses/strains/displacements in the whole domain. It could be the qp or some other points. If so, probably more things we can do in the postprocessing.
Try Term.get_physical_qps() with your term instances, or something like:
from sfepy.discrete.common.mappings import get_physical_qps qps = get_physical_qps(omega, integral) coors = qps.get_merged_values()
The QP coordinates are given in the order of the mesh elements.
For the coordinates of the field nodes (for the lagrange basis), you can use FEField.get_coor() with your field instances (e.g. u.field.get_coor()). For linear fields covering the whole domain (as in your case), the coordinates are the same as the mesh vertices (domain.mesh.coors or mesh.coors).
Does that help?
r.
Yes. that helps a lot. I will test it.
Thanks.
ouyang