SfePy August 2017

sfepy@python.org

3 participants
7 discussions

Calculate total force from cauchy stress

by Patricia Garcia Cañadilla

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

2 years, 6 months

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.

3 years, 4 months

Correct element orientation for tetrahedral mesh

by Kathrin Sobe

Hi Robert, I was able to load a mesh from a nastran .bdf file. The mesh is build of tetrahedrons (nodes: 12593, tetrahedrons: 310106). When I create the Domain it stops with the following error: sfepy: warning: bad element orientation, trying to correct... sfepy: warning: bad element orientation, trying to correct... Traceback (most recent call last): runtime, nodes = tool_displacement() domain = FEDomain('domain', mesh) self.fix_element_orientation() raise RuntimeError('elements cannot be oriented! (%s)' % key) RuntimeError: elements cannot be oriented! (3_4) What is a correct element orientation for tetrahedral elements? The Nastran file was generated automatically. Thanks for your help, Kathrin

3 years, 11 months

Fwd: SfePy - Electro-thermo-mechanics problem

by Robert Cimrman

Hi Manuel, I am forwarding your message to the mailing list, because I could not reach you using your address (Connection timed out error). I will answer to the list as well. r. -------- Forwarded Message -------- Subject: SfePy - Electro-thermo-mechanics problem Date: Mon, 28 Aug 2017 12:31:33 +0200 From: Manuel Feuchter <manuel.feuchter(a)twt-gmbh.de> To: cimrman3(a)ntc.zcu.cz Hey Robert, first of all: SfePy is an amazing project! I am getting to like SfePy more and more :) I tried to start a thread at https://mail.python.org/mm3/archives/list/sfepy@python.org But unfortunately am I not able to. I recieve a "Server Error (500)" . So I try it via Email :) Concerning SfePy: I get stuck with a few problems/questions: Background: I would like to set up a subsequent coupled system (so not fully coupled): 1. electric field evolves an electric current (at the end of the day, I would like to apply the electric field as an time (time step) dependent function.) 2. electric current results in liberated heat by joule heating 3. the liberated heat drives the thermal expansion und results in thermal stresses # 1: ------ Concerning example: http://sfepy.org/doc-devel/examples/multi_physics/thermal_electric.html (multi_physics/thermal_electric.py) Equation: dw_laplace.2.Omega( m.electric_conductivity, psi, phi ) = 0 Your comment states: """ First solve the stationary electric conduction problem. Then use its results to solve the evolutionary heat conduction problem. """ Alright. However, just to understand it comprehensively? When you say you solve the "stationary electric conduction problem". Does this imply actually two equations? $(1) \textbf{E} = -\nabla \phi$\\ $(2) \textbf{j} = \frac{1}{\varrho}\cdot\textbf{E}$ Is this right? # 2: ------ Assuming I got it right, I would like to solve the electric conduction problem in terms of a quasi-static circumstance (external electric potential and internal resistivity will change). This means I would like to solve the electric conduction problem for every time step that I run the simulation and hence use the result (of $\textbf{j}$) for every next time step. I tried so, however I must have gotten sth wrong. The apparent electric field is not used?! Am I right? And how am I able to achieve so? (please find my code below) # 3: ------ For the heat conduction equation: ... = dw_electric_source.2.Omega( m.electric_conductivity, s, phi ) Here, for the source term $\textbf{p} = \varrho \cdot \textbf{j}^2$ is solved? And how does the specific heat gets involved in your equation? Is the source term divided by the thermal conductivity [W/mK]? Does specific_heat to be named exactly like this to be involved in the 'dw_volume_dot.2.Omega( s, dT/dt )' term? # 4: ------ How am I able to apply a time (time step) dependent function for the electric potential?! Looking forward for your reply! Cheers and kind regards from Germany Manuel -------------------------------------------------------- def stress_strain_ext(out, pb, state, extend=False): """ Calculate and output strain and stress for given displacements. """ from sfepy.base.base import Struct from sfepy.mechanics import tensors ev = pb.evaluate strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg') stress = ev('ev_cauchy_stress.2.Omega(m.mech_solid, u)', mode='el_avg') vms = tensors.get_von_mises_stress(stress.squeeze()) vms.shape = (vms.shape[0], 1, 1, 1) out['cauchy_strain'] = Struct(name='output_data', mode='cell', data=strain, dofs=None) out['cauchy_stress'] = Struct(name='output_data', mode='cell', data=stress, dofs=None) out['von_mises_stress'] = Struct(name='output_data', mode='cell', data=vms, dofs=None) return out #! Presettings #! ------- t0 = 0.0 t1 = 100.0 n_step = 11 # material parameters and constants lam = 10.0 mu = 5.0 thermal_expandability = 1.25e-5 T0 = 20.0 # Background temperature. specific_heat = 1.2 thermal_kappa_xyz = 2.0 electric_sigma_xyz = 1.5 #! Options #! ------- options = { 'ts': 'ts', 'nls': 'newton', 'ls': 'ls', 'output_dir' : output_dir, 'post_process_hook' : 'stress_strain_ext', } #! Regions #! ------- regions = { 'Omega': 'all', 'Omega_Cell_Surface': ('vertices of surface', 'facet'), 'Mech_Mount': ('vertices in (x < 0.05)', 'facet'), 'Mech_Front': ('vertices in (x > 0.95)', 'facet'), 'Therm_inlet': ('vertices in (x < 90.05) & (x > 85.95) & (y > 9.95)', 'facet'), 'Therm_outlet': ('vertices in (x < 0.05)', 'facet'), 'Elec_Phi_0' : ('vertices in (x < 37.05) & (x > 32.95) & (y > 9.95)', 'facet'), 'Elec_Phi_1' : ('vertices in (x < 90.05) & (x > 85.95) & (y < 0.05)', 'facet'), } #! Materials #! --------- eye_sym = np.array([[1], [1], [1], [0], [0], [0]], dtype=np.float64) materials = { 'm': ({ 'thermal_conductivity': thermal_kappa_xyz, 'electric_conductivity': electric_sigma_xyz, 'mech_solid' : stiffness_from_lame(3, lam=lam, mu=mu), 'alpha_therm' : (3.0 * lam + 2.0 * mu) * thermal_expandability * eye_sym },), 'mech_load': ({'.val': [0.0, 0.0, 0.0]},), } #! Fields #! ------ fields = { 'displacement': ('real', 'vector', 'Omega', 1), 'temperature': ('real', 1, 'Omega', 1), 'potential' : ('real', 1, 'Omega', 1), } #! Variables #! --------- variables = { 'u': ('unknown field', 'displacement', 2), 'v': ('test field', 'displacement', 'u'), 'T': ('unknown field', 'temperature', 0, 1), # 1 means, set history "on" 's': ('test field', 'temperature', 'T'), 'phi': ('unknown field', 'potential', 1, 1), 'psi': ('test field', 'potential', 'phi'), } #! Initial Conditions #! ------------------- ics = { 'ic' : ('Omega', {'u.all' : 0.0, 'T.0': 0.0, 'phi.all': 0.0, }), } #! Boundary Conditions #! ------------------- ebcs = { 'u0': ('Mech_Mount', {'u.all': 0.0}), 't0': ('Mech_Mount', {'T.0': 0.0}), 't1': ('Mech_Front', {'T.0': 0.0}), 'p0' : ('Elec_Phi_0', {'phi.0' : 0.0}), 'p1' : ('Elec_Phi_1', {'phi.0' : 100.0}), } #! Equations #! --------- equations = { 'potential' : """dw_laplace.2.Omega( m.electric_conductivity, psi, phi ) = 0""", 'heat_conduct' : """dw_volume_dot.2.Omega( s, dT/dt ) + dw_laplace.2.Omega( m.thermal_conductivity, s, T ) = dw_electric_source.2.Omega( m.electric_conductivity, s, phi )""", 'thermo_mech': """dw_lin_elastic.2.Omega(m.mech_solid, v, u) - dw_biot.2.Omega(m.alpha_therm, v, T) = 0""", } #! Solvers #! ------- solvers = { 'ls' : ('ls.scipy_direct', {}), 'newton' : ('nls.newton', { 'i_max' : 1, 'eps_a' : 1e-10, 'eps_r' : 1.0, 'problem': 'nonlinear', }), 'ts': ('ts.simple', { 't0': t0, 't1': t1, 'dt': None, 'n_step': n_step, }), } -------------------------------------------------------- Schöne Grüße kind regards i.A. Manuel Feuchter Competence Unit Engineering ------------------------------------------------------------------------------------------------------------------------------------- Dr.-Ing. Manuel Feuchter Head of Multibody & Structure Analysis TWT GmbH - Science & Innovation Ernsthaldenstraße 17 D-70565 Stuttgart E-Mail: manuel.feuchter(a)twt-gmbh.de Mobil: +49 162 / 212 48 16 www.twt-gmbh.de ------------------------------------------------------------------------------------------------------------------------------------- Geschäftsführung: Dr. Dimitrios Vartziotis, Joachim Laicher (Stv.), Frank Beutenmüller (Stv.) Registergericht: Amtsgericht Stuttgart, HRB Nr. 212778 Umsatzsteuer: ID-Nr.: DE147841145 -------------------------------------------------------------------------------------------------------------------------------------

3 years, 11 months

Fwd: Re: Re: Several questions for solving a linear elastic problem

by Robert Cimrman

Forwarding to the list... -------- Forwarded Message -------- Subject: Re: [sfepy] Re: Several questions for solving a linear elastic problem Date: Tue, 29 Aug 2017 10:59:10 +0200 From: Kathrin Sobe <kathrinsbriefkasten(a)googlemail.com> To: Robert Cimrman <cimrman3(a)ntc.zcu.cz> Hi Robert, just to let you know: This is the result of the example code after I changed the things according to your suggestions. Thank you. I think the questions are answered and I also get a nice result picture. In the next phase I will try the same with the real mesh object. My modifications are as follows: * the surface region is extracted directly with: surface = domain.create_region('Toppart', 'vertices of surface', 'facet') # 528 facets * the material function returns a force vector for each quadrature point: I just splitted the triangle force vector by the number of quadrature points per triangle (facet). I hope this is the right way to distributet the forces. * result extraction with post processing: this works after adding copy_materials = False: The second approach from examples/linear_elasticity/modal_analysis.py, around line 291 did not work in my example. The code is there in comment lines. The code was working, but the result was zero for all elements? However with post_processing I am able to access the results. * Furthermore, I changed the definition of the linear elastic material to solid = Material('solid', values={'D': stiffness_from_youngpoisson(dim=3, young=200, poisson=0.3)}). In my case the material is defined by the youngs-modulus and the poissions ratio. * the gravity is changed to: t3 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v) The values for the material are also defined by a function. Here is the complete modified code example: # coding: utf-8 from __future__ import print_function from __future__ import absolute_import from argparse import ArgumentParser import dill import functools import numpy import time import sys from sfepy.discrete.common.mappings import get_physical_qps sys.path.append('.') from sfepy import data_dir from sfepy.base.base import IndexedStruct, Struct from sfepy.discrete import (FieldVariable, Material, Integral, Integrals, Equation, Equations, Problem, Function) from sfepy.discrete.fem import Mesh, FEDomain, Field from sfepy.terms import Term from sfepy.discrete.conditions import Conditions, EssentialBC from sfepy.solvers.ls import ScipyDirect from sfepy.solvers.nls import Newton from sfepy.postprocess.viewer import Viewer from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson """u aus: http://sfepy.org/doc-devel/tutorial.html """ help = { 'show' : 'show the results figure', } r""" Linear elasticity with pressure traction load on a surface and constrained to one-dimensional motion. Siehe Beispiel als Skript hier!!!!!!!!!1 http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tracti… Find :math:`\ul{u}` such that: .. math:: \int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u}) = - \int_{\Gamma_{right}} \ul{v} \cdot \ull{\sigma} \cdot \ul{n} \;, \quad \forall \ul{v} \;, where .. math:: D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) + \lambda \ \delta_{ij} \delta_{kl} \;. and :math:`\ull{\sigma} \cdot \ul{n} = \bar{p} \ull{I} \cdot \ul{n}` with given traction pressure :math:`\bar{p}`. """ def post_process(out, problem, state, extend=False): """ This will be called after the problem is solved. Parameters ---------- out : dict The output dictionary, where this function will store additional data. problem : Problem instance The current Problem instance. state : State instance The computed state, containing FE coefficients of all the unknown variables. extend : bool The flag indicating whether to extend the output data to the whole domain. It can be ignored if the problem is solved on the whole domain already. Returns ------- out : dict The updated output dictionary. """ # Cauchy strain averaged in elements. strain = problem.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg') out['cauchy_strain'] = Struct(name='output_data', mode='cell', data=strain, dofs=None) # Cauchy stress averaged in elements. problem.equations.materials[0] stress = problem.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)', mode='el_avg', copy_materials=False) out['cauchy_stress'] = Struct(name='output_data', mode='cell', data=stress, dofs=None) # # size 1782 mal 6 Werte Pro Element eine Zeile (=1782) return out def _get_loadcase(): u""" Returns a list of facets and a list of force vectors. The order defines the matching between facet and force vector. """ facets = dill.load(open(fname, "r")) # #force vectors from another example: temporary solution: take the necessary number of vectors force_vectors_example = dill.load(open(fname2, "r")) force_vectors = force_vectors_example[:len(facets)] return facets, force_vectors def _set_not_top_region_zero(facets, force_vectors, facet_indices_top): u""" Surface facets that are not in the top region get the zero force vector. """ zero_vector = numpy.array([0.0, 0.0, 0.0]) for i, f in enumerate(facets): if f not in facet_indices_top: force_vectors[i] = zero_vector else: pass return facets, force_vectors def _get_qp_force_vectors_from_loadcase(force_vectors, facets, region, qps, number_qps_facet): u""" Ensure order of region_facets and divide force vectors for quadrature points. Number of quadrature points per facets is always the same? Even if the triangles have a different size? """ force_vectors_qps = numpy.repeat(force_vectors, 3, axis=0) force_vectors_qps = force_vectors_qps / float(number_qps_facet) return force_vectors_qps def linear_tension(ts, coors, mode=None, force_vectors_qp=None, **kwargs): u""" mapping between force vectors (one for every triangle) and quadrature points. coors are quadrature points: here are three quadrature points per tetrahedron that means: the force vector is divided by three and set to every qp what happens when the number of quadrature points is not even for all triangles? """ if mode == 'qp': # shape: 128x3 # Array with tension values for each quadrature point (constant -0.8) # val = numpy.tile(-10.0, (coor.shape[0], 1, 1)) # numpy.tile(1.0, (3, 1, 1)): [[[1]], [[1]], [[1]]] # val = numpy.tile([[0.001], [0.01], [1.8]], (coors.shape[0], 1, 1)) # 1584 mal [[0.5 -0.45284157 0.4803068], ...] val = numpy.reshape(force_vectors_qp, (coors.shape[0], 3, 1)) return {'val' : val} def get_gravity(ts, coors, mode=None, **kwargs): u""" Gravity in z-direction. Gravity is a volume load: here constant for every quadrature point. Mapping also necessary (see above) """ if mode == 'qp': val = numpy.tile([[0.0], [0.0], [10.5]], (coors.shape[0], 1, 1)) return {'val' : val} def cube_displacement_medium(): u""" 3D Box with tetrahedron elements. Problem type: linear elastic problem Input: force vectors for the top triangles (top surface of the tetrahedron mesh) and gravity for each tetrahedron. Output: - displacement for each vertex - stress tensor """ parser = ArgumentParser() parser.add_argument('--version', action='version', version='%(prog)s') parser.add_argument('-s', '--show', action="store_true", dest='show', default=True, help=help['show']) options = parser.parse_args() # load mesh mesh = Mesh.from_file(data_dir + '/meshes/3d/cube_medium_tetra.mesh') print("Mesh number of nodes: {}, number of tetrahedrons: {}".format(mesh.n_nod, mesh.n_el)) domain = FEDomain('domain', mesh) # why refine? # domain.refine() omega = domain.create_region('Omega', 'all') # Gesamtbereich # facets: bounding box here: -0.5, -0.5, -0.5 bis 0.5, 0.5, 0.5 bb = mesh.get_bounding_box() top = domain.create_region('Top', 'vertices in (z >= 0.48)', 'facet') # 88 facets # gives you vertices of each face in the top region, that can be used to index mesh.coors (mesh.coors[vof]). vertices_indexes_facets_top = mesh.cmesh.get_incident(0, top.faces, 2).reshape((-1, 3)) # vertices_facets_top = mesh.coors[vertices_indexes_facets_top] # facet_indices_top = top.get_facet_indices() surface = domain.create_region('Toppart', 'vertices of surface', 'facet') # 528 facets # vertices of facets # defines the order of facets: order has to be applied on facets load vectors vertices_indexes_facets_surface = mesh.cmesh.get_incident(0, surface.faces, 2).reshape((-1, 3)) # gives you vertices of each face in the top region, that can be used to index mesh.coors (mesh.coors[vof]). z_min = bb[0][2] z_max = bb[1][2] z_range = z_max - z_min z_tol = 0.01 * z_range z_min_tol = z_min + z_tol bottom = domain.create_region('Bottom', 'vertices in (z <= {})'.format(z_min_tol), 'vertex') field = Field.from_args('displacement', numpy.float64, 3, omega, approx_order=1) u = FieldVariable('u', 'unknown', field, order=0) v = FieldVariable('v', 'test', field, primary_var_name='u') # solid = Material('solid', values={'D': stiffness_from_lame(dim=3, lam=5.769, mu=3.846)}) # changed to the given matarial parameter e-modulus and poisson solid = Material('solid', values={'D': stiffness_from_youngpoisson(dim=3, young=200, poisson=0.3)}) # young=29 oder 200 # # loadcase with forces on the surface # # currently: 528 elements and 528 force vectors facets, force_vectors = _get_loadcase() facets, force_vectors = _set_not_top_region_zero(facets, force_vectors, vertices_indexes_facets_top) integral = Integral('i', order=2) # # quadrature points in the order of surface facets # # number of quadrature points for each facet defines the splitting factor? # # e.g. three qps per facet: divide force vector by three # # returns a force vector array with one force vector for each qp qps1 = get_physical_qps(surface, integral) qps_shape = qps1.shape # (528, 3, 3) number_qps_facet = qps_shape[1] force_vectors_qp = _get_qp_force_vectors_from_loadcase(force_vectors, facets, surface, qps1, number_qps_facet) force_vectors_qp = force_vectors_qp * numpy.array([0, 0, -1]) linear_tension_function = functools.partial(linear_tension, force_vectors_qp=force_vectors_qp) load = Material('load', function=linear_tension_function) # material defined by values or function # volume forces f (gravity modelling) f = Material('f', function=Function('get_gravity', get_gravity)) t1 = Term.new('dw_lin_elastic(solid.D, v, u)', integral, omega, solid=solid, v=v, u=u) t2 = Term.new('dw_surface_ltr(load.val, v)', integral, surface, load=load, v=v) # old: # t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid, v=v, u=u) t3 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v) # ???? is this right?? eq1 = Equation('balance_of_forces', t1 + t2 + t3) eqs = Equations([eq1]) fixbottom = EssentialBC('fixbottom', bottom, {'u.all' : 0.0}) ls = ScipyDirect({}) nls_status = IndexedStruct() # configurable residuum # absolute_tolerance = 1e-10 absolute_tolerance = 0.1 relative_tolerance = 1.0 nls = Newton({ 'i_max' : 1, 'eps_a' : absolute_tolerance, 'eps_r' : relative_tolerance, 'macheps' : 1e-16, # Linear system error < (eps_a * lin_red). 'lin_red' : 1e-2, 'ls_red' : 0.1, 'ls_red_warp' : 0.001, 'ls_on' : 1.1, 'ls_min' : 1e-5, 'check' : 0, 'delta' : 1e-6, }, lin_solver=ls, status=nls_status) pb = Problem('linear_elastic_traction', equations=eqs, nls=nls, ls=ls) # pb.set_materials([load, solid]) pb.save_regions_as_groups('regions') # ebcs: elastic boundary conditions pb.time_update(ebcs=Conditions([fixbottom])) start_time = time.time() vec = pb.solve() run_time = time.time()-start_time print("runtime medium problem {} {}: {}".format(mesh.n_nod, mesh.n_el, run_time)) print(nls_status) # von mises stress berechnen: is working now!!! copy_materials = False in function! # pb.save_state('linear_elastic_traction_cylinder.vtk', vec, post_process_hook=post_process) name = 'linear_elastic_traction_medium_box.vtk' pb.save_state(name, vec, post_process_hook=post_process) ######### stress and strain results are zero ####### # # alternative for calculating the stress tensor # state = pb.create_state() # aux = state.create_output_dict() # strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', # mode='el_avg', verbose=False) # # # aux['cauchy_strain'] = Struct(name='output_data', # mode='cell', data=strain, # dofs=None) # # # stress = pb.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)', # mode='el_avg', verbose=False, copy_materials=False) # # aux['cauchy_stress'] = Struct(name='output_data', # mode='cell', data=stress, # dofs=None) # # size 1782 mal 6 Werte Pro Element eine Zeile (=1782) if options.show: view = Viewer(name) view(vector_mode='warp_norm', rel_scaling=2, is_scalar_bar=True, is_wireframe=True) pb.reset() u.reset() v.reset() return run_time, mesh.n_nod if __name__ == '__main__': runtime_small_cube, nodes_small_cube = cube_displacement_medium() 2017-08-25 17:03 GMT+02:00 Robert Cimrman <cimrman3(a)ntc.zcu.cz>: > Hi Kathrin, > > > On 08/25/2017 04:42 PM, Kathrin Sobe wrote: > >> Hi Robert, >> >> thanks a lot for your help. >> >> >> 2017-08-25 14:37 GMT+02:00 Robert Cimrman <cimrman3(a)ntc.zcu.cz>: >> >> Hi Kathrin, >>> >>> On 08/25/2017 10:37 AM, Kathrin Sobe wrote: >>> >>> Hi Robert, >>>> >>>> I try to solve a linear elastic problem. I give a example in a python >>>> file that is a modified version of the example given in >>>> http://sfepy.org/doc-devel/examples/linear_elasticity/linear >>>> _elastic_tractions.html >>>> >>>> The problem I want to solve is the following: >>>> I have a object, defined as a 3D mesh of tetrahedrons. In the test >>>> example I take the medium size 3D cube. The real 3D object is much >>>> bigger >>>> and has a more structured surface. >>>> >>>> On the top of the object there are certain forces. These forces are >>>> defined as 3D vectors, one for each surface triangle. Additionally the >>>> object has certain material parameters given by density, poisson ratio >>>> and >>>> elastic modulus. >>>> >>>> Working through the linear elastic examples I figured out two different >>>> types of solutions: >>>> 1. modeling the forces at the surface as a material with a function >>>> (manual 4.3.9) >>>> 2. modeling the forces at the surface as a boundary condition >>>> >>>> So the first question is in general: Is this the right way to model my >>>> problem? At the moment I think solution type number one is the better >>>> approach? >>>> >>>> >>> Yes. The dw_surface_ltr term you use corresponds to the Neumann boundary >>> conditions = applied forces/tractions/pressure in the weak form of the >>> elasticity PDE. Because you need a position-dependent value of the load, >>> you need to use a material parameter given by a function. >>> >> >> >> ok, than I will concentrate on that approach. Additionally, I want to >> model >> the gravity. For that I added the rho parameter to the material defintion >> and defined a mass equation like this: >> # # rho ist die Dichte: bei geringerer Dichte, wird der Block mehr >> verformt: ist klar >> solid = Material('solid', values={'D': stiffness_from_lame(dim=3, >> lam=5.769, mu=3.846), 'rho': 0.501}) >> load = Material('load', function=linear_tension) # material defined >> by >> values or function >> >> integral = Integral('i', order=2) >> t1 = Term.new('dw_lin_elastic(solid.D, v, u)', integral, omega, >> solid=solid, v=v, u=u) >> t2 = Term.new('dw_surface_ltr(load.val, v)', integral, surface, >> load=load, v=v) >> t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid, >> v=v, u=u) >> eq1 = Equation('elasticity', t1 - t2) >> eq2 = Equation('mass', t3) >> eqs = Equations([eq1, eq2]) >> >> I am not sure if this is the right way to add the gravity to the problem. >> > > Use the dw_volume_lvf term for volume forces just as you use > dw_surface_ltr for surface forces. The material parameter here is a force > per unit volume. The mass matrix would be useful in a dynamic simulation, > not here. > > >>> Going through the first approach I have the following more detailed >>> >>>> problems: >>>> * Region definition: I need a way to define the surface triangles of the >>>> 3D object. Because of the structure of the object, this cannot be done >>>> with >>>> coordinate ranges. >>>> >>>> >>> Is it going to be applied on the whole surface? Anyway, you can select >>> the >>> whole surface of a mesh with 'vertices of surface', for example: >>> >>> surf = domain.create_region('surf', 'vertices of surface', 'facet') >>> >>> Check the manual section 4.3.3 (on regions) for additional possibilities >>> of defining regions and combining existing regions using various >>> operations. >>> >>> Or, if your mesh comes from some mesh generator/cad software, you can >>> mark >>> various vertex / cell groups there and use those for region definitions. >>> This depends on the mesher/format you use - let me know if you need help >>> with that. >>> >>> >> The surface extraction is exactly what I need, cool. I found that out as >> well. Thanks for clarifying the indexing. I will try that. >> >> Just as an information: Not all surface triangles have a force vector, >> only >> the triangles in the upper region of the geometry. >> > > Then you may need to remove the extra faces from the surface region (or > simply define the force zero there). > > >>> * Material definition with force vectors: the 3D force vectors are given >>> >>>> for each surface triangle (defined at the centroid of each triangle). In >>>> the material function examples I found only the definition of forces at >>>> quadrature points and only as scalar values. >>>> >>>> >>> You can use vectors there too, for example: >>> >>> val = nm.tile([[-0.8], [0.1], [0.2]], (coor.shape[0], 1, 1)) >>> >>> To have a constant value per element, you need to know in which element >>> or >>> face each quadrature point is. This information is available in the >>> material function (for example, kwargs['term'].region.faces contains the >>> list of faces of the region, etc.), but the actual code depends on how >>> your >>> material function is represented - is it a table of (cell-value) pairs, >>> or >>> something else? >>> >>> >> At the moment my force defintion is exacly as you expected: cell-value >> pairs. The cell is the surface triangle and the value is the force vector. >> >> Is it that easy: Find the cell for each quadrature point and return the >> cell force vector as its load? Do I need to split the force vector of the >> cell between its quadrature points? Is there always exacly one cell for >> each quadrature point? >> > > You need to provide a force per unit area vector (so yes, we may call it > splitting) in each quadrature point. Each quadrature point is in exactly > one face, yes. The number of quadrature point per face can be found as > follows: > > qps = t2.get_physical_qps() > print qps.shape > (88, 3, 3) > > = 88 faces, 3 QPs per face, 3 coordinates (in 3D) > > So this tells you, that the coor argument of linear_tension has 3 rows for > each face (in the order of top.faces (= t2.region.faces). > > r. > > _______________________________________________ > SfePy mailing list > sfepy(a)python.org > https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/ >

3 years, 11 months

Re: Several questions for solving a linear elastic problem

by Robert Cimrman

Hi Kathrin, On 08/25/2017 04:42 PM, Kathrin Sobe wrote: > Hi Robert, > > thanks a lot for your help. > > > 2017-08-25 14:37 GMT+02:00 Robert Cimrman <cimrman3(a)ntc.zcu.cz>: > >> Hi Kathrin, >> >> On 08/25/2017 10:37 AM, Kathrin Sobe wrote: >> >>> Hi Robert, >>> >>> I try to solve a linear elastic problem. I give a example in a python >>> file that is a modified version of the example given in >>> http://sfepy.org/doc-devel/examples/linear_elasticity/linear >>> _elastic_tractions.html >>> >>> The problem I want to solve is the following: >>> I have a object, defined as a 3D mesh of tetrahedrons. In the test >>> example I take the medium size 3D cube. The real 3D object is much bigger >>> and has a more structured surface. >>> >>> On the top of the object there are certain forces. These forces are >>> defined as 3D vectors, one for each surface triangle. Additionally the >>> object has certain material parameters given by density, poisson ratio and >>> elastic modulus. >>> >>> Working through the linear elastic examples I figured out two different >>> types of solutions: >>> 1. modeling the forces at the surface as a material with a function >>> (manual 4.3.9) >>> 2. modeling the forces at the surface as a boundary condition >>> >>> So the first question is in general: Is this the right way to model my >>> problem? At the moment I think solution type number one is the better >>> approach? >>> >> >> Yes. The dw_surface_ltr term you use corresponds to the Neumann boundary >> conditions = applied forces/tractions/pressure in the weak form of the >> elasticity PDE. Because you need a position-dependent value of the load, >> you need to use a material parameter given by a function. > > > ok, than I will concentrate on that approach. Additionally, I want to model > the gravity. For that I added the rho parameter to the material defintion > and defined a mass equation like this: > # # rho ist die Dichte: bei geringerer Dichte, wird der Block mehr > verformt: ist klar > solid = Material('solid', values={'D': stiffness_from_lame(dim=3, > lam=5.769, mu=3.846), 'rho': 0.501}) > load = Material('load', function=linear_tension) # material defined by > values or function > > integral = Integral('i', order=2) > t1 = Term.new('dw_lin_elastic(solid.D, v, u)', integral, omega, > solid=solid, v=v, u=u) > t2 = Term.new('dw_surface_ltr(load.val, v)', integral, surface, > load=load, v=v) > t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid, > v=v, u=u) > eq1 = Equation('elasticity', t1 - t2) > eq2 = Equation('mass', t3) > eqs = Equations([eq1, eq2]) > > I am not sure if this is the right way to add the gravity to the problem. Use the dw_volume_lvf term for volume forces just as you use dw_surface_ltr for surface forces. The material parameter here is a force per unit volume. The mass matrix would be useful in a dynamic simulation, not here. >> >> Going through the first approach I have the following more detailed >>> problems: >>> * Region definition: I need a way to define the surface triangles of the >>> 3D object. Because of the structure of the object, this cannot be done with >>> coordinate ranges. >>> >> >> Is it going to be applied on the whole surface? Anyway, you can select the >> whole surface of a mesh with 'vertices of surface', for example: >> >> surf = domain.create_region('surf', 'vertices of surface', 'facet') >> >> Check the manual section 4.3.3 (on regions) for additional possibilities >> of defining regions and combining existing regions using various operations. >> >> Or, if your mesh comes from some mesh generator/cad software, you can mark >> various vertex / cell groups there and use those for region definitions. >> This depends on the mesher/format you use - let me know if you need help >> with that. >> > > The surface extraction is exactly what I need, cool. I found that out as > well. Thanks for clarifying the indexing. I will try that. > > Just as an information: Not all surface triangles have a force vector, only > the triangles in the upper region of the geometry. Then you may need to remove the extra faces from the surface region (or simply define the force zero there). >> >> * Material definition with force vectors: the 3D force vectors are given >>> for each surface triangle (defined at the centroid of each triangle). In >>> the material function examples I found only the definition of forces at >>> quadrature points and only as scalar values. >>> >> >> You can use vectors there too, for example: >> >> val = nm.tile([[-0.8], [0.1], [0.2]], (coor.shape[0], 1, 1)) >> >> To have a constant value per element, you need to know in which element or >> face each quadrature point is. This information is available in the >> material function (for example, kwargs['term'].region.faces contains the >> list of faces of the region, etc.), but the actual code depends on how your >> material function is represented - is it a table of (cell-value) pairs, or >> something else? >> > > At the moment my force defintion is exacly as you expected: cell-value > pairs. The cell is the surface triangle and the value is the force vector. > > Is it that easy: Find the cell for each quadrature point and return the > cell force vector as its load? Do I need to split the force vector of the > cell between its quadrature points? Is there always exacly one cell for > each quadrature point? You need to provide a force per unit area vector (so yes, we may call it splitting) in each quadrature point. Each quadrature point is in exactly one face, yes. The number of quadrature point per face can be found as follows: qps = t2.get_physical_qps() print qps.shape (88, 3, 3) = 88 faces, 3 QPs per face, 3 coordinates (in 3D) So this tells you, that the coor argument of linear_tension has 3 rows for each face (in the order of top.faces (= t2.region.faces). r.

3 years, 11 months

Several questions for solving a linear elastic problem

by Kathrin Sobe

Hi Robert, I try to solve a linear elastic problem. I give a example in a python file that is a modified version of the example given in http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tracti… The problem I want to solve is the following: I have a object, defined as a 3D mesh of tetrahedrons. In the test example I take the medium size 3D cube. The real 3D object is much bigger and has a more structured surface. On the top of the object there are certain forces. These forces are defined as 3D vectors, one for each surface triangle. Additionally the object has certain material parameters given by density, poisson ratio and elastic modulus. Working through the linear elastic examples I figured out two different types of solutions: 1. modeling the forces at the surface as a material with a function (manual 4.3.9) 2. modeling the forces at the surface as a boundary condition So the first question is in general: Is this the right way to model my problem? At the moment I think solution type number one is the better approach? Going through the first approach I have the following more detailed problems: * Region definition: I need a way to define the surface triangles of the 3D object. Because of the structure of the object, this cannot be done with coordinate ranges. * Material definition with force vectors: the 3D force vectors are given for each surface triangle (defined at the centroid of each triangle). In the material function examples I found only the definition of forces at quadrature points and only as scalar values. * After solving the problem I need to extract the node displacement and the stress. For the stress calculation I tried to use a post_process function as visible in the example. But this was not working. Thank you for your help. I hope the code in my example code makes the problems more clear. Regards, Kathrin ******************** code example: solution type with force vectors defined in material parameter ******************* # coding: utf-8 from __future__ import print_function from __future__ import absolute_import from argparse import ArgumentParser import numpy as nm import time import sys sys.path.append('.') from sfepy.base.base import IndexedStruct from sfepy.discrete import (FieldVariable, Material, Integral, Equation, Equations, Problem, Function) from sfepy.discrete.fem import Mesh, FEDomain, Field from sfepy.terms import Term from sfepy.discrete.conditions import Conditions, EssentialBC from sfepy.solvers.ls import ScipyDirect from sfepy.solvers.nls import Newton from sfepy.postprocess.viewer import Viewer from sfepy.mechanics.matcoefs import stiffness_from_lame """u aus: http://sfepy.org/doc-devel/tutorial.html """ help = { 'show' : 'show the results figure', } r""" Linear elasticity with pressure traction load on a surface and constrained to one-dimensional motion. Siehe Beispiel als Skript hier!!!!!!!!!1 http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tracti… Find :math:`\ul{u}` such that: .. math:: \int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u}) = - \int_{\Gamma_{right}} \ul{v} \cdot \ull{\sigma} \cdot \ul{n} \;, \quad \forall \ul{v} \;, where .. math:: D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) + \lambda \ \delta_{ij} \delta_{kl} \;. and :math:`\ull{\sigma} \cdot \ul{n} = \bar{p} \ull{I} \cdot \ul{n}` with given traction pressure :math:`\bar{p}`. """ def post_process(out, problem, state, extend=False): """ This will be called after the problem is solved. Parameters ---------- out : dict The output dictionary, where this function will store additional data. problem : Problem instance The current Problem instance. state : State instance The computed state, containing FE coefficients of all the unknown variables. extend : bool The flag indicating whether to extend the output data to the whole domain. It can be ignored if the problem is solved on the whole domain already. Returns ------- out : dict The updated output dictionary. """ from sfepy.base.base import Struct # Cauchy strain averaged in elements. strain = problem.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg') out['cauchy_strain'] = Struct(name='output_data', mode='cell', data=strain, dofs=None) # Cauchy stress averaged in elements. stress = problem.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)', mode='el_avg') out['cauchy_stress'] = Struct(name='output_data', mode='cell', data=stress, dofs=None) return out # # coor: 264 x3 --> qp points in top region! def linear_tension(ts, coor, mode=None, **kwargs): if mode == 'qp': # shape: 128x3 # Array with tension values for each quadrature point (constant -0.8) val = nm.tile(-0.8, (coor.shape[0], 1, 1)) # numpy.tile(1.0, (3, 1, 1)): [[[1]], [[1]], [[1]]] # ## with vectors it does not work # val = nm.tile(nm.array([-0.01, 0.05, -0.8]), (coor.shape[0], 1, 1)) return {'val' : val} def cube_displacement_medium(): u""" 3D Box with tetrahedron elements. Problem type: linear elastic problem linear element approach (later quadratic) Input: I have force vectors for the top triangles (top surface of the tetrahedron mesh) and gravity for each tetrahedron. Output: - displacement for each vertex - stress tensor """ from sfepy import data_dir parser = ArgumentParser() parser.add_argument('--version', action='version', version='%(prog)s') parser.add_argument('-s', '--show', action="store_true", dest='show', default=True, help=help['show']) options = parser.parse_args() # load mesh mesh = Mesh.from_file(data_dir + '/meshes/3d/cube_medium_tetra.mesh') print("Mesh number of nodes: {}, number of tetrahedrons: {}".format(mesh.n_nod, mesh.n_el)) # Domain mit Mesh definieren: hier 448 Knoten, 1782 Elemente domain = FEDomain('domain', mesh) # why refine? domain.refine() omega = domain.create_region('Omega', 'all') # Gesamtbereich # facets: bounding box here: -0.5, -0.5, -0.5 bis 0.5, 0.5, 0.5 bb = mesh.get_bounding_box() top = domain.create_region('Top', 'vertices in (z >= 0.48)', 'facet') bottom = domain.create_region('Bottom', 'vertices in (z <= -0.48)', 'facet') field = Field.from_args('displacement', nm.float64, 3, omega, approx_order=1) u = FieldVariable('u', 'unknown', field, order=0) v = FieldVariable('v', 'test', field, primary_var_name='u') solid = Material('solid', values={'D': stiffness_from_lame(dim=3, lam=5.769, mu=3.846)}) load = Material('load', function=linear_tension) # material defined by values or function integral = Integral('i', order=2) t1 = Term.new('dw_lin_elastic(solid.D, v, u)', integral, omega, solid=solid, v=v, u=u) t2 = Term.new('dw_surface_ltr(load.val, v)', integral, top, load=load, v=v) eq = Equation('elasticity', t1 - t2) eqs = Equations([eq]) fixbottom = EssentialBC('fixbottom', bottom, {'u.all' : 0.0}) ls = ScipyDirect({}) nls_status = IndexedStruct() nls = Newton({ 'i_max' : 1, 'eps_a' : 1e-10, 'eps_r' : 1.0, 'macheps' : 1e-16, # Linear system error < (eps_a * lin_red). 'lin_red' : 1e-2, 'ls_red' : 0.1, 'ls_red_warp' : 0.001, 'ls_on' : 1.1, 'ls_min' : 1e-5, 'check' : 0, 'delta' : 1e-6, }, lin_solver=ls, status=nls_status) pb = Problem('linear_elastic_traction', equations=eqs, nls=nls, ls=ls) pb.save_regions_as_groups('regions') # ebcs: elastic boundary conditions pb.time_update(ebcs=Conditions([fixbottom])) start_time = time.time() vec = pb.solve() run_time = time.time()-start_time print("runtime medium problem {} {}: {}".format(mesh.n_nod, mesh.n_el, run_time)) print(nls_status) # von mises stress calculation (not working) pb.save_state('linear_elastic_simulation2.vtk', vec) # pb.save_state('linear_elastic_traction_cylinder.vtk', vec, post_process_hook=post_process) if options.show: view = Viewer('linear_elastic_simulation2.vtk') view(vector_mode='warp_norm', rel_scaling=2, is_scalar_bar=True, is_wireframe=True) pb.reset() u.reset() v.reset() return run_time, mesh.n_nod if __name__ == '__main__': runtime_small_cube, nodes_small_cube = cube_displacement_medium()

3 years, 11 months

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

SfePy August 2017