Mailman 3 December 2017 - SfePy

Material Parameters of Individual Nodes
by derj...＠web.de 12 Apr '22

12 Apr '22

Hello! It's a great tool you built there! I'd like to use it for some kind of topology optimisation. Is there an "idiomatic" way to set the material parameters of individual nodes? Thank you very much! JonnyB

3 5

Calculate total force from cauchy stress
by Patricia Garcia Cañadilla 17 Jan '19

17 Jan '19

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

4 5

Re: Torque
by Robert Cimrman 10 Mar '18

10 Mar '18

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.

4 7

ANN: SfePy 2017.4
by Robert Cimrman 29 Dec '17

29 Dec '17

I am pleased to announce release 2017.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 or by the isogeometric analysis (limited support). It is distributed under the new BSD license. Home page: http://sfepy.org Mailing list: https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/ Git (source) repository, issue tracker: https://github.com/sfepy/sfepy Highlights of this release -------------------------- - basic support for penalty-based contacts - support for user-defined contexts in all solvers and preconditioners - new example: dispersion analysis of heterogeneous periodic materials For full release notes see http://docs.sfepy.org/doc/release_notes.html#id1 (rather long and technical). Cheers, Robert Cimrman --- Contributors to this release in alphabetical order: Robert Cimrman Jan Heczko Lubos Kejzlar Jan Kopacka Vladimir Lukes Matyas Novak

1 0

Extracting physical surfaces/volumes named in gmsh
by maciejz4＠hotmail.com 08 Dec '17

08 Dec '17

Hello, Is there a possibility of extracting physical groups from .msh file somehow? I work on some complex mesh and that would help me a lot with setting boundary condition and material constraints.

4 3

Re: Direct solver without recalculating matrix factorization
by Robert Cimrman 06 Dec '17

06 Dec '17

Hi David, the problem was that you re-created the whole problem in each loading step. Check the attached modified script, and try playing with 'presolve' and 'is_linear' settings. IMHO in your case it might be easier to just use a problem description file, where you define a quasi-static time-stepper, and a load by a function - see sample2.py (adapted from [1] to compute the same problem as sample.py). But with the interactive version you can easily add further features. r. [1] http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tracti… On 12/05/2017 04:53 PM, David Fernández wrote: > I have attached the code sample and the modified ls.py which I used for > making sure that the factorization doesn't repeat. > > Best, > > David > > > On 05/12/17 14:09, Robert Cimrman wrote: >> Hi David, >> >> could you send me the modified example? >> >> r. >> >> On 12/05/2017 02:04 PM, David Fernández wrote: >>> Hi, thank you for your reply, it has been very useful. However, there >>> are still some things which are a bit strange to me. >>> >>> 1- I have modified >>> "http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_intera…" >>> >>> to solve the system many times in a loop, updating only the body forces >>> (see code bellow) and setting ls = >>> ScipyDirect({'method':'umfpack','presolve':True}) outside the loop (and >>> not defining nls solver). The problem is that the factorization still >>> takes place in every iteration. I have looked for the reason and it lies >>> in line 163 of ls.py (if self.mtx_id != id(mtx)). For some reason the >>> system matrix is being updated even if I'm not mollifying it. By >>> manually enforcing to no repeat the factorization (modifying ls.py) it >>> works and it's much faster. Is it a bug or am I doing something wrong in >>> the sample code? >>> >>> 2- Using nls = Newton({'is_linear' : True}, lin_solver=ls, >>> status=nls_status) (commented lines in the sample code bellow) the >>> results are not good (the residuals are huge and the solution are bad) >>> >>> >>> Best, >>> >>> David. >>> >>> >>> #!/usr/bin/env python >>> from __future__ import print_function >>> from __future__ import absolute_import >>> from argparse import ArgumentParser >>> import numpy as nm >>> >>> import sys >>> sys.path.append('.') >>> >>> from sfepy.base.base import IndexedStruct >>> from sfepy.discrete import (FieldVariable, Material, Integral, Function, >>> Equation, Equations, Problem) >>> 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 >>> import time >>> >>> def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0): >>> """ >>> Define a displacement depending on the y coordinate. >>> """ >>> val = shift * coors[:,1]**2 >>> >>> return val >>> >>> helps = { >>> 'show' : 'show the results figure', >>> } >>> >>> def main(): >>> 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=False, help=helps['show']) >>> options = parser.parse_args() >>> >>> mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh') >>> domain = FEDomain('domain', mesh) >>> >>> min_x, max_x = domain.get_mesh_bounding_box()[:,0] >>> eps = 1e-8 * (max_x - min_x) >>> omega = domain.create_region('Omega', 'all') >>> gamma1 = domain.create_region('Gamma1', >>> 'vertices in x < %.10f' % (min_x + >>> eps), >>> 'facet') >>> gamma2 = domain.create_region('Gamma2', >>> 'vertices in x > %.10f' % (max_x - >>> eps), >>> 'facet') >>> >>> field = Field.from_args('fu', nm.float64, 'vector', omega, >>> approx_order=4) >>> >>> u = FieldVariable('u', 'unknown', field) >>> v = FieldVariable('v', 'test', field, primary_var_name='u') >>> >>> m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0)) >>> >>> fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0}) >>> >>> bc_fun = Function('shift_u_fun', shift_u_fun, >>> extra_args={'shift' : 0.01}) >>> shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun}) >>> >>> ls = ScipyDirect({'method':'umfpack','presolve':True}) >>> t = time.time() >>> for i in range(100): >>> >>> f = Material('f', val=[[0.02+i*0.001], [0.01]]) >>> integral = Integral('i', order=4) >>> t1 = Term.new('dw_lin_elastic(m.D, v, u)', >>> integral, omega, m=m, v=v, u=u) >>> t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, >>> f=f, v=v) >>> eq = Equation('balance', t1 + t2) >>> eqs = Equations([eq]) >>> #nls_status = IndexedStruct() >>> #nls = Newton({'is_linear' : True}, lin_solver=ls, >>> status=nls_status) >>> #pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls) >>> >>> pb = Problem('elasticity', equations=eqs, ls=ls) >>> >>> pb.time_update(ebcs=Conditions([fix_u, shift_u])) >>> vec = pb.solve() >>> pb.save_state('/home/me/linear_elasticity_{}.vtk'.format(i), >>> vec) >>> >>> print(time.time()-t) >>> if __name__ == '__main__': >>> main() >>> >>> >>> >>> On 04/12/17 22:22, Robert Cimrman wrote: >>>> Hello David, >>>> >>>> On 12/04/2017 03:01 PM, David Fernández wrote: >>>>> Hello, I'm very new to sfepy and therefore it's being a bit hard to >>>>> find >>>>> in the documentation a solution to my problem, which first of all I >>>>> would like to know if it is solvable within sfepy. >>>>> >>>>> I'm solving the elastic equilibrium equations, which I need to solve >>>>> many times during a single program run. However, in between calls >>>>> to the >>>>> solver, only the RHS of the FEM formulation is changing, i.e. the body >>>>> forces or Neumann conditions. As in this case the stiffness matrix is >>>>> constant, I would like to use a direct solver (e.g. umfpack) to >>>>> calculate the matrix factorization (i.e. setup the solver) only >>>>> once, an >>>>> then simply call the solver many times with different body forces. I >>>>> would really appreciate if some of you could point me in the right >>>>> direction, or provide me with a couple of lines of code showing how >>>>> this >>>>> can be (more or less) achieved with this really nice module. >>>> >>>> Yes, this could be done basically in two ways: either modify one of >>>> the non-stationary examples with a time-stepper (e.g. [1]) where you >>>> set the option >>>> >>>> 'quasistatic' : True, >>>> >>>> for the 'ts.simple' time-stepper, and the option >>>> >>>> 'is_linear' : True, >>>> >>>> for the 'nls.newton' nonlinear solver and the option >>>> >>>> 'presolve' : True, >>>> >>>> for the 'ls.scipy_direct' linear solver. >>>> this will result in prefactorization of the matrix before the first >>>> time step, and reuse of the factors in all the time steps. >>>> >>>> The other way is to take a script such as [2] and modify it as needed. >>>> >>>> Feel free to ask more. >>>> Cheers, >>>> r. >>>> >>>> [1] >>>> http://sfepy.org/doc-devel/examples/diffusion/time_advection_diffusion.html >>>> >>>> [2] >>>> http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_intera… >>>> >>>> _______________________________________________ >>>> SfePy mailing list >>>> sfepy(a)python.org >>>> https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/ >>> >> >

1 0

Re: Direct solver without recalculating matrix factorization
by Robert Cimrman 05 Dec '17

05 Dec '17

Hi David, could you send me the modified example? r. On 12/05/2017 02:04 PM, David Fernández wrote: > Hi, thank you for your reply, it has been very useful. However, there > are still some things which are a bit strange to me. > > 1- I have modified > "http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_intera…" > to solve the system many times in a loop, updating only the body forces > (see code bellow) and setting ls = > ScipyDirect({'method':'umfpack','presolve':True}) outside the loop (and > not defining nls solver). The problem is that the factorization still > takes place in every iteration. I have looked for the reason and it lies > in line 163 of ls.py (if self.mtx_id != id(mtx)). For some reason the > system matrix is being updated even if I'm not mollifying it. By > manually enforcing to no repeat the factorization (modifying ls.py) it > works and it's much faster. Is it a bug or am I doing something wrong in > the sample code? > > 2- Using nls = Newton({'is_linear' : True}, lin_solver=ls, > status=nls_status) (commented lines in the sample code bellow) the > results are not good (the residuals are huge and the solution are bad) > > > Best, > > David. > > > #!/usr/bin/env python > from __future__ import print_function > from __future__ import absolute_import > from argparse import ArgumentParser > import numpy as nm > > import sys > sys.path.append('.') > > from sfepy.base.base import IndexedStruct > from sfepy.discrete import (FieldVariable, Material, Integral, Function, > Equation, Equations, Problem) > 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 > import time > > def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0): > """ > Define a displacement depending on the y coordinate. > """ > val = shift * coors[:,1]**2 > > return val > > helps = { > 'show' : 'show the results figure', > } > > def main(): > 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=False, help=helps['show']) > options = parser.parse_args() > > mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh') > domain = FEDomain('domain', mesh) > > min_x, max_x = domain.get_mesh_bounding_box()[:,0] > eps = 1e-8 * (max_x - min_x) > omega = domain.create_region('Omega', 'all') > gamma1 = domain.create_region('Gamma1', > 'vertices in x < %.10f' % (min_x + eps), > 'facet') > gamma2 = domain.create_region('Gamma2', > 'vertices in x > %.10f' % (max_x - eps), > 'facet') > > field = Field.from_args('fu', nm.float64, 'vector', omega, > approx_order=4) > > u = FieldVariable('u', 'unknown', field) > v = FieldVariable('v', 'test', field, primary_var_name='u') > > m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0)) > > fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0}) > > bc_fun = Function('shift_u_fun', shift_u_fun, > extra_args={'shift' : 0.01}) > shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun}) > > ls = ScipyDirect({'method':'umfpack','presolve':True}) > > t = time.time() > for i in range(100): > > f = Material('f', val=[[0.02+i*0.001], [0.01]]) > > integral = Integral('i', order=4) > > t1 = Term.new('dw_lin_elastic(m.D, v, u)', > integral, omega, m=m, v=v, u=u) > t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v) > eq = Equation('balance', t1 + t2) > eqs = Equations([eq]) > > #nls_status = IndexedStruct() > #nls = Newton({'is_linear' : True}, lin_solver=ls, > status=nls_status) > #pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls) > > pb = Problem('elasticity', equations=eqs, ls=ls) > > pb.time_update(ebcs=Conditions([fix_u, shift_u])) > vec = pb.solve() > pb.save_state('/home/me/linear_elasticity_{}.vtk'.format(i), vec) > > print(time.time()-t) > > if __name__ == '__main__': > main() > > > > On 04/12/17 22:22, Robert Cimrman wrote: >> Hello David, >> >> On 12/04/2017 03:01 PM, David Fernández wrote: >>> Hello, I'm very new to sfepy and therefore it's being a bit hard to find >>> in the documentation a solution to my problem, which first of all I >>> would like to know if it is solvable within sfepy. >>> >>> I'm solving the elastic equilibrium equations, which I need to solve >>> many times during a single program run. However, in between calls to the >>> solver, only the RHS of the FEM formulation is changing, i.e. the body >>> forces or Neumann conditions. As in this case the stiffness matrix is >>> constant, I would like to use a direct solver (e.g. umfpack) to >>> calculate the matrix factorization (i.e. setup the solver) only once, an >>> then simply call the solver many times with different body forces. I >>> would really appreciate if some of you could point me in the right >>> direction, or provide me with a couple of lines of code showing how this >>> can be (more or less) achieved with this really nice module. >> >> Yes, this could be done basically in two ways: either modify one of >> the non-stationary examples with a time-stepper (e.g. [1]) where you >> set the option >> >> 'quasistatic' : True, >> >> for the 'ts.simple' time-stepper, and the option >> >> 'is_linear' : True, >> >> for the 'nls.newton' nonlinear solver and the option >> >> 'presolve' : True, >> >> for the 'ls.scipy_direct' linear solver. >> this will result in prefactorization of the matrix before the first >> time step, and reuse of the factors in all the time steps. >> >> The other way is to take a script such as [2] and modify it as needed. >> >> Feel free to ask more. >> Cheers, >> r. >> >> [1] >> http://sfepy.org/doc-devel/examples/diffusion/time_advection_diffusion.html >> [2] >> http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_intera… >> _______________________________________________ >> SfePy mailing list >> sfepy(a)python.org >> https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/ >

1 0

Direct solver without recalculating matrix factorization
by David Fernández 05 Dec '17

05 Dec '17

Hello, I'm very new to sfepy and therefore it's being a bit hard to find in the documentation a solution to my problem, which first of all I would like to know if it is solvable within sfepy. I'm solving the elastic equilibrium equations, which I need to solve many times during a single program run. However, in between calls to the solver, only the RHS of the FEM formulation is changing, i.e. the body forces or Neumann conditions. As in this case the stiffness matrix is constant, I would like to use a direct solver (e.g. umfpack) to calculate the matrix factorization (i.e. setup the solver) only once, an then simply call the solver many times with different body forces. I would really appreciate if some of you could point me in the right direction, or provide me with a couple of lines of code showing how this can be (more or less) achieved with this really nice module. Thank very much!

2 1

Re: Problems with edited Script
by Robert Cimrman 03 Dec '17

03 Dec '17

Hi Dio, On 11/23/2017 03:22 PM, Dio Dwianto wrote: > Hi Robert, > > thanks for your help! That finally solved it, apparently the initial > condition was set to distribute temperature between -2 and 2 degrees while > I set the temperature at 20 degrees Celcius. > > Now I'm trying to add another boundary condition to the equation to take > into account the convective heat transfer that's going on the surface. I > tried to use dw_surface_dot, but I'm not sure which parameter I shoud use. > This is what I wrote so far: > > [image: Inline image 1] > > with coef.val I meant the convective heat transfer coefficient [W/m^2K] The dw_surface_dot needs to be used as follows: equations = { 'Temperature' : """dw_volume_dot.i.Omega(Omega.capacity, s, dT/dt) + dw_laplace.i.Omega(Omega.conductivity, s, T) - dw_surface_dot.i.Omega(coef.val, s, T) = 0""" } Note also, that there is no ds/dt in the dw_volume_dot term. > But all I get was an error message: > > [image: Inline image 2] > > Can you tell me what I should do to fix this? Also do I have to add another > solver to solve the convective heat transfer problem? You might want to use dw_surface_integrate term to apply a given heat flux on a boundary - see [1] for a usage example. No other solver needs to be added. [1] http://sfepy.org/doc-devel/examples/diffusion/poisson_neumann.html r.

1 0