I have some problems by applying a load force by means of a function. I have checked available information, and tested many ways without success. So, please, could you help me?
The problem is the next. The problem has 5 steps, I want to apply the force in a specific nodes, defined as forze_zone, proportionally to the step number. I think my problem is about the output format of the get_load function. Here is the code.
def get_load(ts, coors, mode='qp',equations=None, term=None, problem=None, **kwargs):
if mode == 'qp':
val= nm.empty(coors.shape[0], dtype=nm.float64)
val=14000*ts.step/ts.n_step/n_coors #force nodes array format [index nodes, ]
val.reshape(coors.shape[0], 1, 1) #output format (# of nodes, 1,1)
return {'.val': val}
functions = {
'get_pars' : (get_pars,), # funcion used to get material properties
'get_load' : (get_load,),
'get_zone_fix': (_get_zone_fix,), # function to get nodes to fix the part
'get_zone_force': (_get_zone_force,), #function to get nodes to applied the load
}
options = {
'ts' : 'ts',
'output_format' : 'vtk',
'save_times' : 'all',
'post_process_hook' : 'post_process',
}
mesh = Mesh.from_file(filename_mesh)
domain = FEDomain('domain', mesh)#Create a domain
aux_coors=nm.copy(mesh.coors) # copy of mesh coordinates
n_coors=len(_get_zone_force(aux_coors)) # number of nodes in load zone
force_nodes=_get_zone_force(aux_coors) #node index to apply the force
gdata = geometry_data['3_4']
nc = len(gdata.coors)
integral_vn = Integral('ivn', coors=gdata.coors, weights=[gdata.volume / nc] * nc)
regions = {
'Omega' : 'all',
'Fix_zone' : ('vertices by get_zone_fix', 'vertex'),
'Force_zone' : ('vertices by get_zone_force','vertex'),
}
materials = {
'nonlinear' : 'get_pars',
'Load' : 'get_load',
}
fields = {
'displacement': ('real', 'vector', 'Omega', 1),
}
variables = {
'u' : ('unknown field', 'displacement', 0),
'v' : ('test field', 'displacement', 'u'),
}
ebcs = {
'Fixed' : ('Fix_zone', {'u.all' : 0.0}),
}
integrals = {
'i': 2, # 2nd order quadrature over a 3 dimensional space.
}
equations = {
'balance_of_forces in time' :
"""dw_lin_elastic.2.Omega(nonlinear.D, v, u) =
dw_point_load.1.Force_zone(Load.val, v)""",
}
solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton',
{'i_max' : 1,
'eps_a' : 1e-10,
'eps_r' : 1.0,
}),
'ts' : ('ts.simple',
{'t0' : 0.0,
't1' : 1.0,
'dt' : None,
'n_step' : 5,
'quasistatic' : True,
'verbose' : 1,
}),
}
Thank You
Jesus Sanz
