Mailman 3 June 2021 - 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

ANN: SfePy 2021.2
by Robert Cimrman 15 Jul '21

15 Jul '21

I am pleased to announce the release of SfePy 2021.2. Description ----------- SfePy (simple finite elements in Python) is a software for solving systems of coupled partial differential equations by finite element methods. It is distributed under the new BSD license. Home page: https://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 -------------------------- - new sensitivity analysis terms - positive FE basis based on Bernstein polynomials - smaller memory footprint of terms with constant material parameters For full release notes see [1]. Cheers, Robert Cimrman [1] http://docs.sfepy.org/doc/release_notes.html#id1 --- Contributors to this release in alphabetical order: Robert Cimrman Vladimir Lukes

2 1

Element stiffness matrix
by 979439736＠qq.com 15 Jul '21

15 Jul '21

I would like to ask what is the law of arrangement of the stiffness matrix of the element? I look forward to your reply, thank you!

4 9

Material coefficients with time derivatives of field variables
by Max Maximilian 12 Jul '21

12 Jul '21

Hey all, I found my script cannot successfully compute the volume_dot matrix including time derivatives of field variables, since it stayed 0 at all time. Does it mean I need to reformulate the equations that exclude time derivative terms? ## CODE from __future__ import absolute_import from sfepy.examples.dg.example_dg_common import * from sfepy import data_dir from sfepy.base.base import output import numpy as nm from sfepy.discrete.fem import Mesh from sfepy.discrete.fem.meshio import UserMeshIO import matplotlib.pyplot as plt from sfepy.base.base import Struct filename_mesh=None X1 = 0 XN = 2e-4 n_nod = 2001 n_el = n_nod - 1 if filename_mesh is None: filename_mesh = get_gen_1D_mesh_hook(X1, XN, n_nod) approx_order=2 t0=0.0 t1=0.5 porosity = 0.52 # PTL porosity theta_PTL = nm.pi*30./180. # PTL contact angle KG = 4.5e-15 # Permeability, m2 R = 8.3145 # Ideal gas constant, J/mol/K T = 273.15+80 # Ambient temperature, K patm = 101325.0 # Ambient pressure sCh = 0.5 # Channel liquid saturation sl = 0.4 # Left boundary liquid saturation Mox = 0.032 # Molar mass density of oxygen, kg/mol i_app = 12000. # Applied current density, A/m2 F = 96485 # Faraday constant, C/mol Mw = 0.018 # Molar mass density of water, kg/mol muG = 1.881e-5 # Gas viscosity, Pa*s muL = 3.56e-4 # Liqiud viscosity, Pa*s rho_H2O = 996 # Water density, kg/m3 m = 3 # Relative permeabiltiy coefficient psat = 10**(8.07131-1730.63/(233.426+T-273.15))*133.322 # Saturated water gas pressure, Pa kc = 1e4 # Condensation rate coefficient, 1/s kv = 5.1e-5 # Evaporation rate coefficient, 1/(Pa*s) Dox = 2.1e-5 # Diffusivity of oxygen, m^2/s gas_rate = i_app/4/F*Mox liquid_rate = -i_app/2/F*Mw # Boundary & Initial pressure/mass fraction pc0 = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-sCh)-2.120*(1-sCh)**2+1.263*(1-sCh)**3) rho00 = Mox*(patm+pc0-psat)/R/T+psat*Mw/R/T c0 = Mox*(patm+pc0-psat)/(Mox*(patm+pc0-psat)+psat*Mw) pcl = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-0.4)-2.120*(1-0.4)**2+1.263*(1-0.4)**3) rhol = Mox*(patm+pcl-psat)/R/T+psat*Mw/R/T cl = Mox*(patm+pcl-psat)/(Mox*(patm+pcl-psat)+psat*Mw) example_name = "fluid_dynamics_1D" def get_ic_w(coors, ic=None): pc = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))-2.120*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))**2+1.263*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))**3) rho0 = Mox*(patm+pc-psat)/R/T+psat*Mw/R/T rho0.shape = (coors.shape[0], 1, 1) return rho0 def get_ic_s(coors, ic=None): s0 = (sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl s0.shape = (coors.shape[0], 1, 1) return s0 def get_ic_c(coors, ic=None): pc = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))-2.120*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))**2+1.263*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))**3) c0 = (Mox*(patm+pc-psat))/(Mox*(patm+pc-psat)+psat*Mw) c0.shape = (coors.shape[0], 1, 1) return c0 def post_process(out, pb, state, extend=False): # Field variable evalution rho_grad = pb.evaluate('ev_grad.i.Omega(rho)', mode='el_avg', verbose=False) s_grad = pb.evaluate('ev_grad.i.Omega(s)', mode='el_avg', verbose=False) c_grad = pb.evaluate('ev_grad.i.Omega(c)', mode='el_avg', verbose=False) liquid_saturation_values = pb.evaluate('ev_volume_integrate.i.Omega(s)', mode='el_avg', verbose=False) rho_values = pb.evaluate('ev_volume_integrate.i.Omega(rho)', mode='el_avg', verbose=False) conc_values = pb.evaluate('ev_volume_integrate.i.Omega(c)', mode='el_avg', verbose=False) dpcds = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(-1.417+2.120*2*(1-liquid_saturation_values)-3*1.263*(1-liquid_saturation_values)**2) liquid_velocity = -KG*liquid_saturation_values**m/muL*R*T*((conc_values/Mox+(1-conc_values)/Mw)*rho_grad+(1/Mox-1/Mw)*rho_values*c_grad-dpcds*s_grad/R/T) gas_velocity = -KG*(1-liquid_saturation_values)**m/muG*R*T*((conc_values/Mox+(1-conc_values)/Mw)*rho_grad+(1/Mox-1/Mw)*rho_values*c_grad) out['liquid_velocity'] = Struct(name='output_data', mode='cell', data=liquid_velocity, dofs=None) out['gas_velocity'] = Struct(name='output_data', mode='cell', data=gas_velocity, dofs=None) return out def get_coef(ts, coors, problem, equations=None, mode=None, **kwargs): """ Calculates the coefficient matrix and returns them. This relation results in velocities and diffusivity. """ if mode == 'qp': # Field and gradient values in quadrature points coordinates given by integral i # - they are the same as in `coors` argument. if ts.step == 0: liquid_saturation_values = (sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl else: liquid_saturation_values = problem.evaluate('ev_volume_integrate.i.Omega(s)', mode='qp', verbose=False) liquid_saturation_values.shape = (coors.shape[0], 1, 1) if ts.step == 0: liquid_saturation_dt = nm.zeros(coors.shape[0]) else: liquid_saturation_dt = problem.evaluate('ev_volume_integrate.i.Omega(ds/dt)', mode='qp', verbose=False) liquid_saturation_dt.shape = (coors.shape[0], 1, 1) if ts.step == 0: dpcds = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(-1.417+2.120*2*(1-liquid_saturation_values)-3*1.263*(1-liquid_saturation_values)**2) grad_rho_values = Mox/R/T*dpcds*(sCh-sl)/(XN-X1) else: grad_rho_values = problem.evaluate('ev_grad.i.Omega(rho)', mode='qp', verbose=False) grad_rho_values.shape = (coors.shape[0], 1, 1) if ts.step == 0: pc = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-liquid_saturation_values)-2.120*(1-liquid_saturation_values)**2+1.263*(1-liquid_saturation_values)**3) rho_values = Mox*(patm+pc-psat)/R/T+psat*Mw/R/T else: rho_values = problem.evaluate('ev_volume_integrate.i.Omega(rho)', mode='qp', verbose=False) rho_values.shape = (coors.shape[0], 1, 1) if ts.step == 0: rho_dt = nm.zeros(coors.shape[0]) else: rho_dt = problem.evaluate('ev_volume_integrate.i.Omega(drho/dt)', mode='qp', verbose=False) rho_dt.shape = (coors.shape[0], 1, 1) if ts.step == 0: conc_values = (Mox*(patm+pc-psat))/(Mox*(patm+pc-psat)+psat*Mw) else: conc_values = problem.evaluate('ev_volume_integrate.i.Omega(c)', mode='qp', verbose=False) conc_values.shape = (coors.shape[0], 1, 1) if ts.step == 0: grad_conc_values = psat*Mw/(Mox*(patm+pc-psat)+psat*Mw)**2*Mox*dpcds*(sCh-sl)/(XN-X1) else: grad_conc_values = problem.evaluate('ev_grad.i.Omega(c)', mode='qp', verbose=False) grad_conc_values.shape = (coors.shape[0], 1, 1) # Gas velocity u1 = -KG*(1-liquid_saturation_values)**m/muG*R*T*((conc_values/Mox+(1-conc_values)/Mw)*grad_rho_values+(1/Mox-1/Mw)*rho_values*grad_conc_values) ## Gas density mass conservation # Mass matrix of gas density: (1-s)*eps mass_dens = (1-liquid_saturation_values)*porosity # Source matrix of gas density: porosity*ds/dt source_dens = liquid_saturation_dt*porosity # Diffusion coupling matrix of gas density: u1 diffcop_dens = u1 # Reaction matrix of gas density: Rsat (Phase transistion) p = rho_values*(1-conc_values)*R*T/Mw+conc_values*rho_values*R*T/Mox reac_dens = kc*porosity*(1-liquid_saturation_values)*Mw/R/T*(rho_values*(1-conc_values)*R*T/Mw-psat)*nm.heaviside(rho_values*(1-conc_values)*R*T/Mw-psat, 0) + kv*porosity*liquid_saturation_values*rho_H2O*(rho_values*(1-conc_values)*R*T/Mw-psat)*nm.heaviside(psat-rho_values*(1-conc_values)*R*T/Mw, 0) ## Liquid water mass conservation # Mass matrix of liquid saturation: porosity*rho_H2O mass_sat = rho_H2O*porosity*nm.ones(coors.shape[0]) mass_sat.shape = (coors.shape[0], 1, 1) # Diffusion_r matrix of liquid saturation: rho_H2O*K*s^m/muL*R*T*((c/Mox+(1-c)/Mw)*grad_rho+grad_c*(1/Mox-1/Mw)*rho) Diff_r_sat = rho_H2O*KG*liquid_saturation_values**m/muL*R*T*((conc_values/Mox+(1-conc_values)/Mw)*grad_rho_values+(1/Mox-1/Mw)*grad_conc_values*rho_values) # Stiffness matrix of liquid saturation: rho_H2O*dpcds*K*s^m/muL dpcds = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(-1.417+2.120*2*(1-liquid_saturation_values)-3*1.263*(1-liquid_saturation_values)**2) stiff_sat = rho_H2O*KG*liquid_saturation_values**m/muL*dpcds # Reaction matrix of liquid saturation: Rsat reac_sat = reac_dens ## Mass contration of Oxygen # Mass matrix of oxygen concentration mass_conc = porosity*(1-liquid_saturation_values)*rho_values # Source 1 matrix of oxygen concentration: eps*(1-s)*drho/dt source1_conc = rho_dt*porosity*(1-liquid_saturation_values) # Source 2 matrix of oxygen concentration: eps*ds/dt*rho source2_conc = liquid_saturation_dt*porosity*rho_values # Diffusion coupling matrix of oxygen concentration: u1*rho diffcop_conc = u1*rho_values # Stiffness matrix of oxygen concentration: rho*Deff stiff_conc = rho_values*Dox*(porosity*(1-liquid_saturation_values))**1.5 output('min:', source_dens.min(), 'max:', source_dens.max()) return {'mass_dens' : mass_dens, 'source_dens' : source_dens, 'diffcop_dens' : diffcop_dens, 'reac_dens' : reac_dens, 'mass_sat' : mass_sat, 'Diff_r_sat' : Diff_r_sat, 'stiff_sat' : stiff_sat, 'reac_sat' : reac_sat, 'mass_conc' : mass_conc, 'source1_conc' : source1_conc, 'source2_conc' : source2_conc, 'diffcop_conc' : diffcop_conc, 'stiff_conc' : stiff_conc} materials = { 'coef' : 'get_coef', 'flux' : ({ 'gas_rate' : gas_rate, 'liquid_rate' : liquid_rate, 'conc_rate' : 0, },), } fields = { 'gas_density': ('real', 'scalar', 'Omega', approx_order), 'liquid_saturation': ('real', 'scalar', 'Omega', approx_order), 'gas_mass_fraction': ('real', 'scalar', 'Omega', approx_order), } variables = { 'rho': ('unknown field', 'gas_density', 0, 1), 'v1': ('test field', 'gas_density', 'rho'), 's': ('unknown field', 'liquid_saturation', 1, 1), 'v2': ('test field', 'liquid_saturation', 's'), 'c': ('unknown field', 'gas_mass_fraction', 2, 1), 'v3': ('test field', 'gas_mass_fraction', 'c'), } regions = { 'Omega' : 'all', 'Gamma_Left' : ('vertices in (x < 1e-7)', 'facet'), 'Gamma_Right' : ('vertices in (x > 1.999e-4)', 'facet'), } ebcs = { 'right' : ('Gamma_Right', {'rho.0' : rho00, 's.0' : sCh, 'c.0' : c0}), 'left' : ('Gamma_Left', {'rho.0' : rhol, 's.0' : 0.4, 'c.0' : cl}), } integrals = { 'i' : 2*approx_order, } equations = { 'gas_density' : """dw_volume_dot.i.Omega(coef.mass_dens, v1, drho/dt) - dw_volume_dot.i.Omega(coef.source_dens, v1, rho) - dw_diffusion_coupling.i.Omega(coef.diffcop_dens, v1, rho) + dw_volume_integrate.i.Omega(coef.reac_dens, v1) = 0""", 'liquid_saturation' : """dw_volume_dot.i.Omega(coef.mass_sat, v2, ds/dt) + dw_diffusion_r.i.Omega(coef.Diff_r_sat, v2) - dw_laplace.i.Omega(coef.stiff_sat, v2, s) - dw_volume_integrate.i.Omega(coef.reac_sat, v2) = 0""", 'gas_mass_fraction' : """dw_volume_dot.i.Omega(coef.mass_conc, v3, dc/dt) + dw_volume_dot.i.Omega(coef.source1_conc, v3, c) - dw_volume_dot.i.Omega(coef.source2_conc, v3, c) - dw_diffusion_coupling.i.Omega(coef.diffcop_conc, v3, c) + dw_laplace.i.Omega(coef.stiff_conc, v3, c) = 0 """, } # -dw_surface_ndot.i.Gamma_Left(flux.gas_rate, v1) -dw_surface_ndot.i.Gamma_Left(flux.liquid_rate, v2) -dw_surface_ndot.i.Gamma_Left(flux.conc_rate, v3) functions = { 'get_coef' : (get_coef,), 'get_ic_w' : (get_ic_w,), 'get_ic_s' : (get_ic_s,), 'get_ic_c' : (get_ic_c,), } ics = { 'ic' : ('Omega', {'rho.0' : 'get_ic_w', 's.0' : 'get_ic_s', 'c.0' : 'get_ic_c'}), } solvers = { 'ls' : ('ls.mumps', { }), 'newton' : ('nls.newton', { 'eps_a' : 1e-4, 'eps_r' : 1e-3, 'i_max' : 5, }), 'tss' : ('ts.adaptive', { 't0' : t0, 't1' : t1, 'dt' : None, 'verbose' : 1, 'quasistatic' : False, 'n_step' : 100, }), } options = { 'ts' : 'tss', 'nls' : 'newton', 'ls' : 'ls', 'save_times' : 10, 'output_dir' : 'output/' + example_name, 'output_format' : "vtk", 'post_process_hook' : 'post_process', }

3 6

Implementation of IGA surface terms
by hugo.levy＠gmail.com 28 Jun '21

28 Jun '21

Hi Robert, I have been using Sfepy for the past two months now, and the more I use it the more I like it! Thank you for developing this great FEM (& more) software. I spent some time trying the various functionalities sfepy offers and recently discovered the IGA mode, which I managed to use interactively. After some testing, it turns out isogeometric analysis might be well suited for my applications, where geometry plays a crucial role. So far, my IGA simulations outperform FEM ones, which really sparked my interest for this type of analysis. Having no background in IGA, I read references [1] and [2] to grasp the mathematical concepts at stake, and reference [3] for the Bézier extraction method. That being said, user's guide mentions that surface terms in IGA are not implemented yet in Sfepy, and I was wondering why? For instance, I would like to compute integrals over (D-1)-dimensional spaces (where D denotes the problem spatial dimension), but as expected, running the line: >>> pb.evaluate('ev_surface_integrate.3.Gamma(u)', mode='eval') results with >>> ValueError: unknown dof connectivity type! (surface) This error is raised in sfepy\discrete\iga\fields.py function setup_extra_data. Of course, I could get around this missing functionality by sampling the scalar field of interest over the desired surface/curve manually (i.e. with field.evaluate_at) and use scipy.integrate. But that would mean giving up on the 'exact' geometry description, which I am not keen on. Is there a particular reason for not implementing those surface terms in IGA, like a technical difficulty that I am not yet aware of? Otherwise, could you give me insights into where to start if I were to implement my own quadrature rule for computing surface integrals in IGA mode? I haven't figured out in which module the evaluation of Bernstein basis functions at Gauss quadrature points occurs. I know juggling with NURBS, Bézier meshes, etc. doesn't seem particularly easy. Yet I'm willing to give it a try. Thank you very much for your help. Best regards, Hugo L. ------------------ [1] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, Elsevier, 2005, 194 (39-41), pp.4135-4195. 10.1016/j.cma.2004.10.008. hal-01513346 [2] Cimrman, Robert. (2014). Enhancing SfePy with Isogeometric Analysis. [3] Michael J. Borden, Michael A. Scott, John A. Evans, Thomas J. R. Hughes: Isogeometric finite element data structures based on Bezier extraction of NURBS, Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, Texas, March 2010.

2 10