On 06/09/2016 09:02 AM, Patricia Garcia Cañadilla wrote:

Hi Robent,

Yes, it was one oh the error, thanks. Now, my problem is I do not how to perform the T = stress \cdot n. and to finally obtain the force (sorry I am very new with python). I have this:

ev = problem.evaluate field = problem.fields['displacement'] region = problem.domain.regions['Gamma'] integrals = problem.get_integrals()['i2'] n = get_normals(field,integrals,region)

Try using:

geo, _ = field.get_mapping(region, integral, 'surface')

then geo contains both normals and the jacobians, and you can integrate it just by something like

np.sum(np.dot(stress, geo.normal) * geo.det)

but first you need to restore the stress to matrix form - it is stored as a vector with 11, 22, 12 components (in 2D). For that you can use [2].

Also you might want to convert the stress to Cauchy stress first [1] - the 2nd Piola-Kirchoff stress that is used in the total Lagrangian formulation has no direct physical meaning.

And there is a catch, see below.

[1] http://sfepy.org/doc-devel/src/sfepy/mechanics/tensors.html?highlight=cauchy... [2] http://sfepy.org/doc-devel/src/sfepy/mechanics/tensors.html?highlight=get_fu...

stress_aux = 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'); T = stress_aux*n; Force = ev('ev_surface_integrate.2.Gamma(T)')

and the shapes are stress_aux = (9699,1,3,1) n=(299,2,2,1)

This is because you evaluate the stress in Omega, while the normals on Gamma. Evaluation of stress on surface is not directly supported for the moment. But one could use some workaround(s) here - could you send me the input file (+ a small mesh, if needed) to try something?

r.

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'),

}

Thanks again :)

El miércoles, 8 de junio de 2016, 17:12:09 (UTC+2), Robert Cimrman escribió:

...

Hi Patricia

...

Dear Robert,

I have a 2D body of hyperelastic material which contracts and I would

...

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

On 06/06/2016 03:27 PM, Patricia Garcia Cañadilla wrote: like 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)

Here you probably meant using 'region' instead of 'regions', right? regions are the dict in the module name space of the problem description file, so the error is pretty easily overlooked :)

If that does not help, try sending a complete file that reproduces the problem.

r.

Hi Andrew,

the reaction forces are the forces in the FE nodes, that are fixed by the Dirichlet boundary conditions - is that what you want?

If yes, such forces can be easily computed by assembling, after the problem has been solved, the residual vector of your equation without the Dirichlet BCs [1], e.g. by:

[1] https://scicomp.stackexchange.com/a/29350

# state is the State containing your variables with the problem solution. state.apply_ebc() nls = problem.get_nls() residual = nls.fun(state())

then the reaction forces in the nodes of a given region can be obtained by:

reg = problem.domain.regions[region_name] dofs = field.get_dofs_in_region(reg, merge=True)

res = residual.reshape((-1, dim)) # dim = space dimension = DOFs per node. reactions = res[dofs]

Some more answers below...

On 1/16/19 5:44 PM, andrew.glusiec@gmail.com wrote:

Not sure if I should open a new thread, but I'm having a very similar problem where I can't figure out how to calculate the reaction forces and torques through a face of a 3D model. I have the model working to the point where I have the following commands execute successfully:

svar = FieldVariable('sigma', 'parameter', field, primary_var_name='sigma') stress = ev('ev_cauchy_stress.3.Omega(m.D, u)', mode='qp', copy_materials=False) svar.set_from_qp(stress, integral) supportstresses = ev('ev_surface_integrate.3.Gamma1(sigma)', var_dict={'sigma':svar}, mode='qp') geo, _ = field.get_mapping(region, integral, 'surface')

I now have what I believe are variables that represent: Cauchy Stress throughout the model (stress)

Yes, in the quadrature points in the volume of the domain. Not surface.

Normal vectors/Jacobians for the face (geo)

Yes.

I am confused as to what supportstresses represents, if it is even a useful quantity, I thought it would be the reaction forces, but the results don't seem to correlate well with the forces I should be seeing

If the field of svar is a surface field, then supportstresses should be the stress interpolated to the surface quadrature points from the surface field nodes of svar. Also, another averaging/interpolation occurred to get the stress from the volume QP to the svar field nodes in svar.set_from_qp(stress, integral).

So I am now lost on how to pick out the stresses that correspond to the face where the support is located and then convert that information to a simple np array for use in this equation:

np.sum(np.dot(stress, geo.normal) * geo.det)

Maybe you are just missing using the region, where you have the support defined? Is it Gamma1?

Also note, that ev_cauchy_strain_s term can be used to directly evaluate the strain on a surface region. The same term for the stress is not implemented, but it would be an easy addition. Or you can just take the strain, and multiply it by the material tensor D_ijkl.

My goal once the reaction forces are working is to also calculate reaction torques if possible based on the distances of each force to some normal axis, but I don't know how I would go about getting the position information needed for this.

I can post my full code if that would be helpful. Does anyone have suggestions that could help determine the reaction forces and moments?

Yes, if the above does not help you, post a minimalistic (small mesh etc.) running example that demonstrates what you want to achieve.

Cheers, r.