Hi Joshua,

I am responding in the text/code below.

On 06/07/2021 17:36, Joshua colletta wrote:

Code Follows:

from argparse import ArgumentParser, RawDescriptionHelpFormatter import numpy as nm import sys import os sys.path.append('.') from sfepy.base.base import IndexedStruct from sfepy.discrete import (FieldVariable, Material, Integral, Equation, Equations, Problem, Functions) 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 from sfepy.base.conf import transform_functions

helps = { 'young': "the Young's modulus [default: %(default)s]", 'poisson': "the Poisson's ratio [default: %(default)s]", 'load': "the vertical load value (negative means compression)" " [default: %(default)s]", 'show': 'show the results figure', 'probe' : 'probe the results', 'order' : 'displacement field approximation order [default: %(default)s]', }

def inner_region_(coors, domain=None): magnet = global_magnet x, y, z = coors[:, 0], coors[:, 1], coors[:, 2] mag_min_x, mag_max_x = (magnet[1][0] - (magnet[0][0] / 2))/ 2, (magnet[1][0] + (magnet[0][0] / 2))/ 2 mag_min_y, mag_max_y = (magnet[1][1] - (magnet[0][1] / 2))/ 2, (magnet[1][1] + (magnet[0][1] / 2))/ 2 mag_min_z, mag_max_z = (magnet[1][2] - (magnet[0][2] / 2))/ 2, (magnet[1][2] + (magnet[0][2] / 2))/ 2 flag = nm.where((mag_max_x > x > mag_min_x) & (mag_max_y > y > mag_min_y) & (mag_max_z > z > mag_min_z))[0]

NumPy arrays cannot be compared using a < b < c syntax, you have to use (a < b) & (b < c):

flag = nm.where((mag_max_x > x) & (x > mag_min_x) & (mag_max_y > y) & (y > mag_min_y) & (mag_max_z > z) & (z > mag_min_z))[0]

 return flag

def simulation(mesh_loc, magnet): parser = ArgumentParser(description=__doc__, formatter_class=RawDescriptionHelpFormatter) parser.add_argument('--version', action='version', version='%(prog)s') parser.add_argument('--load', metavar='float', type=float, action='store', dest='load', default=50.0, help=helps['load']) parser.add_argument('--poisson', metavar='float', type=float, action='store', dest='poisson', default=0.453, help=helps['poisson']) parser.add_argument('--young', metavar='float', type=float, action='store', dest='young', default=470.00, help=helps['young']) parser.add_argument('-s', '--show', action="store_true", dest='show', default=True, help=helps['show']) options = parser.parse_args()

 global global_magnet
 global_magnet = magnet

 conf_functions = {
     'inner_region': (inner_region_,),
 }
 functions = Functions.from_conf(transform_functions(conf_functions))

 mesh = Mesh.from_file(mesh_loc)

 domain = FEDomain('domain', mesh)

 min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
 eps = 0.01 * (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')
 magnet_region = domain.create_region('magnet0', 'vertices by inner_region', 'facet', parent='Omega', functions=functions)

 field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=1)

 u = FieldVariable('u', 'unknown', field)
 v = FieldVariable('v', 'test', field, primary_var_name='u')

 integral = Integral('i', order=3)
 ma = Material('ma', D=stiffness_from_youngpoisson(3, options.young, options.poisson))
 m = Material('m', values={'rho': {'Omega': 2700, 'magnet0': 6000}, 'E': {'Omega': 70e3, 'magnet0': 1.6e5}})

This was just an example of how to define region-dependent values. In your particular case, you can use something like:

 m = Material('m', values={
     'D': {
         'Omega': stiffness_from_youngpoisson(3, options.young, 

options.poisson), 'magnet0': stiffness_from_youngpoisson(3, options.young2, options.poisson2) } })

• the young2, poisson2 options would need to be added.

 load = Material('load', values={'val': options.load})

 t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)
 t3 = Term.new('dw_surface_ltr( load.val, v )', integral, gamma2, load=load, v=v)
 eq = Equation('balance', t1 + t3)

 eqs = Equations([eq])

 fix = EssentialBC('fix', gamma1, {'u.all': 0.0})
 ls = ScipyDirect({})

 nls_status = IndexedStruct()
 nls = Newton({}, lin_solver=ls, status=nls_status)

 pb = Problem('elasticity', equations=eqs)

 pb.save_regions_as_groups('regions')
 pb.set_bcs(ebcs=Conditions([fix]))
 pb.set_solver(nls)
 status = IndexedStruct()
 s = pb.solve(status=status)

 print('Nonlinear solver status:\n', nls_status)
 print('Stationary solver status:\n', status)

 pb.save_state('linear_elasticity.vtk', s)

 if options.show:
     view = Viewer('linear_elasticity.vtk')
     view(vector_mode='warp_norm', rel_scaling=2,
          is_scalar_bar=True, is_wireframe=True)

if __name__ == "__main__": magnet = [[1,1,1],[9,9,9]] ur_path_to_mesh = "" simulation(ur_path_to_mesh, magnet)

I cant seem to attach a mesh, but it takes a 10x10x10 cube as a .msh file!

No problem, I have generated a test mesh using:

python ./script/blockgen.py -d 10,10,10 -c 5,5,5 -s 31,31,31 -o mesh.vtk

nm.where() is still not behaving like i expected. I tried passing it the x,y,z vectors but still no dice.

See the comment in the code above.

Testing the material parameters with the simple nm.where(x<3)[0] lead to a different error, stating that there is no key 'D', I believe stemming from the equations expected use of parameters in:

 t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)

I would like to define my materials using the stiffness_from_youngpoisson() function, is this still possible?

Yes, see above again.

r.