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] 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}})
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!
nm.where() is still not behaving like i expected. I tried passing it the x,y,z vectors but still no dice. 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?