Mailman 3 November 2021 - SfePy

Material Parameters of Individual Nodes
by derj...＠web.de April 11, 2022

April 11, 2022

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

Replicating FEM Study
by Alexander Rockhill Nov. 13, 2021

Nov. 13, 2021

-----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 3

Using Neumann boundary conditions instead of Dirichlet boundary conditions
by contact＠dstoupis.com Nov. 12, 2021

Nov. 12, 2021

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!

3 10