SfePy

sfepy@python.org

852 discussions

ANN: SfePy 2021.4

by Robert Cimrman

I am pleased to announce the release of SfePy 2021.4. 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 -------------------------- - improved pyvista-based visualization script resview.py - gallery images generated using resview.py - homogenization tools: new parallel recovery of multiple microstructures - new "dry water" flow example 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 Yves Delley Florian Le Bourdais Vladimir Lukes Brylie Christopher Oxley Tomáš Zítka

2 weeks, 4 days

Replicating FEM Study

by Alexander Rockhill

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA256 Hello, I'm attempting to replicate the study from Figure 14 here: https://doi.org/10.1016/j.media.2019.06.004 with a regular mesh, not using the meshless method. I have: - A masked brain that can be used to generate a mesh from voxels - The location of points to set for the displacement relative to points on the brain mesh - Coordinates and triangles connecting the coordinates for the inner skull surface I would like to get out of the model: - The displaced vertices so that the original brain image can be resampled What I'm stuck on: - How to define the boundary conditions per coordinate so that the correct coordinates are displaced - How to get the displacement vectors out of the solved problem Bonus: - If there is a way to do this with the inner skull surface acting as a fixed boundary condition frictionless with the brain, that would be preferable, I'm a more lost on how to do this though If someone had time to point me in the right direction a bit, that would be greatly appreciated. I think if I knew a few more of the correct functions to use, I would be very close. Below is what I have so far (requires `pip install mne`). ``` import os.path as op import numpy as np import mne import nibabel as nib from sfepy.base.base import IndexedStruct from sfepy.discrete import (FieldVariable, Material, Integral, Equation, Equations, Problem) from sfepy.discrete.fem import FEDomain, Field from sfepy.mesh.mesh_generators import gen_mesh_from_voxels 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 misc_path = mne.datasets.misc.data_path() # get original and projected positions raw = mne.io.read_raw(op.join(misc_path, 'ecog', 'sample_ecog_ieeg.fif')) trans = mne.coreg.estimate_head_mri_t( 'sample_ecog', op.join(misc_path, 'ecog')) montage = raw.get_montage() montage.apply_trans(trans) # transform to MRI coordinates ch_pos_disp = np.array([ pos for pos in montage.get_positions()['ch_pos'].values()]) raw.info = mne.preprocessing.ieeg.project_sensors_onto_brain( raw.info, trans, 'sample_ecog', subjects_dir=op.join(misc_path, 'ecog')) montage = raw.get_montage() montage.apply_trans(trans) # transform to MRI coordinates ch_pos_orig = np.array([ pos for pos in montage.get_positions()['ch_pos'].values()]) # define model dist = 10 # displace all vertices up to 1 cm from points brainmask = nib.load(op.join(misc_path, 'ecog', 'sample_ecog', 'mri', 'brainmask.mgz')) mesh = gen_mesh_from_voxels(brainmask, [1, 1, 1]) domain = FEDomain('domain', mesh) omega = domain.create_region('Omega', 'all') fixed = domain.create_region( 'Fixed', XXX, 'facet') disp = domain.create_region( 'Displaced', XXX, 'facet') field = Field.from_args('fu', np.float64, 'vector', omega, approx_order=2) u = FieldVariable('u', 'unknown', field) v = FieldVariable('v', 'test', field, primary_var_name='u') brain = Material('Brain', D=stiffness_from_youngpoisson( dim=3, young=3000, poisson=0.49), density=1000) # csf = Material('Brain', D=stiffness_from_youngpoisson( # dim=3, young=100, poisson=0.49), density=1000) integral = Integral('i', order=1) term = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=brain, v=v, u=u) eq = Equation('balance', term) eqs = Equations([eq]) fix_u = EssentialBC('fix_u', fixed, {'u.all': 0.0}) disp_u = EssentialBC( 'disp_u', disp, {'u.0': XXX, 'u.1': XXX, 'u.2': XXX}) ls = ScipyDirect({}) nls_status = IndexedStruct() nls = Newton({}, lin_solver=ls, status=nls_status) pb = Problem('deform', equations=eqs) # view setup ''' pb.save_regions_as_groups('regions') view = Viewer('regions.vtk') view() ''' pb.set_bcs(ebcs=Conditions([fix_u, disp_u])) pb.set_solver(nls) status = IndexedStruct() state = pb.solve(status=status) print('Nonlinear solver status:\n', nls_status) print('Stationary solver status:\n', status) pb.save_state('brain_deform.vtk', state) view = Viewer('brain_deform.vtk') view(vector_mode='warp_norm', rel_scaling=1, is_scalar_bar=True, is_wireframe=True) # modify original image T1 = nib.load(op.join(misc_path, 'ecog', 'sample_ecog', 'mri', 'T1.mgz')) # TO DO: get displacement vector from each vertex out, resample image # based on T1 image intensity of original element positions ``` Thanks, Alex Rockhill -----BEGIN PGP SIGNATURE----- Version: BCPG v1.65 iQEcBAABCgAGBQJhjq3nAAoJEOCF7T6SI2rbi08H/j9PIYQxCaT0aDe5+EmAKmaw I+rZTM+gKBoDZkJD+XTUqUdXKALLxqdpNShi5JAy5MTeXS+CGJ9PufoL3DgDoWJE czfbczqGXzZStsKL7qVROEBF7sP2uvWNxdVYTXqUWVxXdn814SBn+pGHxPF7JShU XzTP1D+Bp/mlPL7/3rcw/qpI0LnFBdLKIcnwEERetEfpjKdj1v0/sAOpAJknHhGO yKwU+4M5bt95Q29S1qyzM5ybMFToj8SEdEBa8qPgtdBt67NlYv4bfrpsgt79i/eq +xgbcWV5M6zDeJyhtlm7YTRhc9cCsxTdug0jB5TdDXWFFsy2kS1NNtacnsy7BVE= =Ntng -----END PGP SIGNATURE-----

2 months

Using Neumann boundary conditions instead of Dirichlet boundary conditions

by contact＠dstoupis.com

Hi everyone, I am trying to solve the following equation: https://www.HostMath.com/Show.aspx?Code=%5Csigma%5Cnabla%5Cphi%2B%5Comega%5… I want to set Neumann boundary conditions on the outer cell group and the condition is this: https://www.HostMath.com/Show.aspx?Code=%5Cnabla%5Cphi%5Ccdot%5Cvec%7Bn%7D%… I have read the previous discussions listed below: https://groups.google.com/g/sfepy-devel/c/8PPndgYC3_4 https://groups.google.com/g/sfepy-devel/c/R9SIMoNiRfc https://github.com/sfepy/sfepy/issues/195 https://groups.google.com/g/sfepy-devel/c/9-cf42ZKsog https://sfepy.org/doc/examples/diffusion-poisson_neumann.html?highlight=neu… and also tried the poisson_neumann example, but I can not understand how to implement a purely Neumann condition solution. I do not understand how to insert the boundary condition in the problem. I tried using the LinearCombinationBC in the following way, although I can not grasp how it works: temporary_domain = domain_definition.create_region('boundary_region', 'r.boundary_Gamma', 'vertex', add_to_regions=False) LinearCombinationBC('boundary_condition', [temporary_domain, None], {'potential.all' : None}, None, 'normal_direction') Below is the code logic I use to define the equations. from sfepy.discrete import Integral, Equation, Equations, Problem from sfepy.terms import Term term_arguments_laplace = { 'c_permitivity': material_definition, 'potential_test': virtual_var, 'potential': state_var, } term_arguments_vector_potential = { 'c_permitivity': material_definition, 'magnetic_vector_potential': parameter_var, 'potential_test': virtual_var, } integral = Integral('i1', order=2) overall_volume = domain_definition.create_region('Omega', 'all') equation_term_laplace = Term.new('dw_laplace(c_permitivity.e, potential_test, potential)', integral, overall_volume, **term_arguments_laplace) equation_term_vector_potential = Term.new('dw_stokes(c_permitivity.e, magnetic_vector_potential, potential_test)', integral, overall_volume, **term_arguments_vector_potential) equations_list.append(Equation('equation', equation_term_laplace + equation_term_vector_potential)) pb = Problem('problem', equations = Equations(equations_list)) If I define Dirichlet conditions, I get a solution successfully, although from the physics point of view it is wrong. Also the solution I get using the Dirichlet conditions is not what I would expect. Can you please give me some directions on how to implement the Neumann conditions, if possible, using SfePy? Thank you!

2 months

updated lagrangian

by Peixun Wang

Hi, I'm trying to solve a large deformation problem interactively using updated lagrangian formulation. I have some trouble in updating the mesh. In the code I attached here, I attempted to define new_ulf_iteration in problem.nls_iter_hook (line 149) and nls.iter_hook (line 152). But both methods failed to update the msh. How do you update mesh in interactive mode? Thanks! Regarls, Peixun

3 months

Integral/Field approximation order

by kshargh.ali＠gmail.com

Hi, My question is about order of approximation in field as well as integral. My understanding is that the order of field is the order of the polynomial shape functions in element and order of integral is for the integration in qp. For a 3 point triangular element as well as 4 point quadrilateral element, I think order of field should be 1 and using any other value is "wrong". Is this correct ? Also, I noticed from user manual that the order of integral = 2* order of field. I was curious what is the basis for this assumption ? Best, Ali

3 months, 1 week

question

by hu Lee

Hi I have some questions about sfepy, the details are in the PPT. Best, Lee

3 months, 1 week

Puzzle

by 979439736＠qq.com

Hi With the same element size, Poisson's ratio and elastic modulus, why is the arrangement order of the element stiffness matrix obtained in ANSYS different from that in sfepy? Can you tell me which finite element book you are referring to? Secondly, I want to set different materials for different finite element elements. What functions do I need to call? Best， Lee

3 months, 1 week

Fwd: Puzzle

by hu Lee

Hi I have some questions about sfepy, the details are in the PPT. Best, Lee

3 months, 1 week

Stress and strain convergence with mesh size

by kshargh.ali＠gmail.com

Dear all, I notice something in my simulations and I wanted to share it with you. I appreciate if I could have your thoughts on it: In one simulation, when I uniaxially strectch a plate and change the mesh size, I noticed that the nodal stresses and strains remain in the same range. However, once I added an elliptical hole in its center, it seems that as I decrease the mesh size, stress and strain especially at the sharp edges of the ellipse goes up and do not reach to a plateau. My understanding is that it shouldn't be real and in fact should reach to a steady state. Please let me know what you think. Best, Ali

3 months, 2 weeks

ANN: SfePy 2021.3

by Robert Cimrman

I am pleased to announce the release of SfePy 2021.3. 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 -------------------------- - unified volume and surface integration terms - improved pyvista-based visualization script 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 Hugo Levy Vladimir Lukes

3 months, 2 weeks

Jump to page:

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

SfePy