On 07/31/2014 02:27 PM, Ouyang wrote:
>
>
> On Thursday, July 31, 2014 6:56:46 PM UTC+8, Robert Cimrman wrote:
>>
>> On 07/31/2014 12:12 PM, Ouyang wrote:
>>> Hi Robert,
>>>
>>>     1. how to output principal stresses in sfepy?  I can calculate
>> 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?)

>>

>>>     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.