Hello, sir
I am trying to implement stress diffusion one way copuling, i am trying to implement time depndent boundary conditions, i am getting issues. Can you please help me in correcting issues i am encountering.
from sfepy import data_dir
from sfepy.discrete import FieldVariable, Material, Integral,Problem
from sfepy.discrete.fem import Mesh
from sfepy.discrete.fem.meshio import UserMeshIO
from sfepy.base.base import Struct
from sfepy.discrete import Equation, Equations
from sfepy.discrete.conditions import EssentialBC
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.terms import Term
from sfepy.solvers import Solver
import numpy as np
from sfepy.discrete.fem import Mesh,FEDomain,Field
from sfepy import terms
from sfepy.base.base import Struct
from sfepy.discrete.fem import Mesh,FEDomain,Field
from sfepy.solvers.ts_solvers import SimpleTimeSteppingSolver
from sfepy.base.base import Struct
# In[702]:
import numpy as np
import meshio # version 4.4.6
from sfepy.discrete import Material
from sfepy.discrete.fem import FEDomain, Field
from sfepy.discrete.fem import Mesh # Sfepy I/O utilities
def converter4sfepy(HKnew, ignoreTags=[], exterior=False, Rcut=None):
mesh = Mesh.from_file(HKnew)
cell_dim = mesh.dim
# All the necessary & sufficient info for defining a mesh
data = list(mesh._get_io_data(cell_dim_only=cell_dim)) # to be mod
# Managing vertex groups
ngroups = np.array([None for _ in range(mesh.n_nod)])
reader = meshio.read(HKnew)
for key in reader.cells_dict.keys():
conn = reader.cells_dict[key]
formerTags = list(reader.cell_data_dict.values())[0][key]
for k, tag in enumerate(formerTags):
if tag not in ignoreTags:
for idx in conn[k]:
if ngroups[idx] is None:
ngroups[idx] = int(tag)
ngroups[np.where(ngroups==None)] = 400 # default marker for untagged
data[1] = ngroups.astype(dtype=np.int32)
# Managing cell groups
conns = list(reader.cells_dict.values())[-1] # entities of highest dim
mat_ids = np.max(ngroups[conns], axis=1)
data[3] = [mat_ids.astype(dtype=np.int32)]
# Save the new mesh
newname = "modified_meshHKnew.vtk"
mesh = Mesh.from_data(newname, *data)
mesh.write(newname, None, binary=False)
out = newname
return out
original_file = 'HKnew.vtk'
modified_file = converter4sfepy(original_file)
mesh = Mesh.from_file(modified_file)
# In[703]:
domain = FEDomain('domain',mesh)
# In[704]:
min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
min_y, max_y = domain.get_mesh_bounding_box()[:, 1]
eps = 1e-8 * (max_x - min_x)
omega = domain.create_region('Omega', 'all')
# Selecting the boundary as the nodes tagged 200 (physical group)
group_id = 29
group_id = 30
Gamma1 = domain.create_region('Left','vertices in x < %.10f' % (min_x + eps),'facet')
Gamma2 = domain.create_region('Right','vertices in x > %.10f' % (max_x - eps),'facet')
Region1 = domain.create_region('Region1', 'cells of group 29')
Region2 = domain.create_region('Region2', 'cells of group 30')
import numpy as np
def get_displacement(ts, coors, bc=None, problem=None):
"""
Define the time-dependent displacement using a sine function.
"""
frequency = 0.2 # Adjust the frequency as needed
amplitude = 0.4 # Adjust the amplitude as needed
time = ts.time
displacement = amplitude * np.sin(2 * np.pi * frequency * time) * coors[:, 0]
return displacement
# In[705]:
def main(cli_args):
dims = parse_argument_list(cli_args.dims, float)
shape = parse_argument_list(cli_args.shape, int)
centre = parse_argument_list(cli_args.centre, float)
material_parameters = parse_argument_list(cli_args.material_parameters,
float)
order = cli_args.order
ts_vals = cli_args.ts.split(',')
ts = {
't0' : float(ts_vals[0]), 't1' : float(ts_vals[1]),
'n_step' : int(ts_vals[2])}
do_plot = cli_args.plot
# In[706]:
import numpy as np
def get_concentration(ts, coors, bc=None, problem=None):
"""
Define the time-dependent concentration using a cosine function.
"""
c_min = 0.01 # Adjust as needed
c_max = 0.99 # Adjust as needed
time = ts.time
concentration = [-(c_max - c_min)/2 * np.cos(np.pi*time/200) + (c_max + c_min)/2]* coors[:, 1]
return concentration
# In[707]:
field_c = Field.from_args('concentration', np.float64, 'scalar', omega , approx_order=2)
field_u = Field.from_args('displacement', np.float64, 'vector', omega,approx_order=1)
# In[708]:
from sfepy.discrete import (FieldVariable, Material, Integral, Function,Equation, Equations, Problem)
c = FieldVariable('c', 'unknown', field_c, history=1)
c_test = FieldVariable('c_test', 'test', field_c, primary_var_name='c')
u = FieldVariable('u', 'unknown', field_u)
v = FieldVariable('v', 'test', field_u, primary_var_name='u')
# In[709]:
integral = Integral('i', order=4)
# In[710]:
from sfepy.discrete.conditions import Conditions, EssentialBC, InitialCondition
# Initial conditions.
def get_ic(coors, ic):
x, y, z = coors.c
return ic_fun
ic_fun = Function('ic_fun', get_ic)
ic = InitialCondition('ic', omega, {'c.0' : 0.01})
# In[711]:
ics = {
'ic' : ('Omega', {'c.0' : 0.01}),
}
# In[720]:
from sfepy.mechanics.matcoefs import stiffness_from_lame
#Steel
lam1 = 65.9340e9 # Lame parameters
mu1 = 76.9320e9
omega1 = 4.17e-6 # Thermal expansion coefficient
#Aluminium
lam2 =2.4360e11 # Lame parameters
mu2 = 2.6717e11
omega2 = 4.17e-6
c0 = 0.0
#This script uses SI units (meters, kilograms, Joules. . . ) except for temperature, which is expressed in degrees Celsius.
eye_sym = np.array([[1], [1], [0]],dtype=np.float64)
m1 = Material('m1', lam1=lam1, mu1=mu1, alpha1 = (3.0 * lam1 + 2.0 * mu1)* omega1 * eye_sym)
m2 = Material('m2', lam2=lam2, mu2=mu2, alpha2 = (3.0 * lam2 + 2.0 * mu2)* omega2 * eye_sym)
# In[713]:
from sfepy import terms
term1 = Term.new('dw_laplace(c_test, c)', integral, Region1,c_test=c_test, c=c)
term2 = Term.new('dw_laplace(c_test, c)', integral, Region2,c_test=c_test, c=c)
term3 = Term.new('dw_lin_elastic_iso(m1.lam1, m1.mu1, v, u)',integral, Region1, m1=m1, v=v, u=u)
term4 = Term.new('dw_lin_elastic_iso(m2.lam2, m2.mu2, v, u)',integral, Region2, m2=m2, v=v, u=u)
term5 = Term.new('dw_biot(m1.alpha1, v, c)',integral, Region1, m1=m1, v=v, c=c)
term6 = Term.new('dw_biot(m2.alpha2, v, c)',integral, Region2, m2=m2, v=v, c=c)
term7 = Term.new('dw_dot(c_test, dc/dt)',integral, Region1, c_test=c_test, c=c)
term8 = Term.new('dw_dot(c_test, dc/dt)',integral, Region2, c_test=c_test, c=c)
eq = Equation('balance', term1*5e-1+ term2*5e-1+ term3 + term4 - term5 - term6 + term7 + term8)
eqs = Equations([eq])
# In[714]:
from sfepy.base.base import IndexedStruct
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
ls = ScipyDirect({})
nls_status = IndexedStruct()
nls = Newton({}, lin_solver=ls, status=nls_status)
# In[715]:
pb = Problem('displacement_concentration', equations=eqs)
# In[716]:
pb.save_regions_as_groups('regions')
# In[717]:
# Assuming Gamma1 and Gamma2 are instances of the Region class
u_left = EssentialBC('u_left', Gamma1, {'u.0': 0.0})
conc_fun = Function('conc_fun', get_concentration)
concentration = EssentialBC('concentration', Gamma2, {'c.0': conc_fun})
ebcs = Conditions([u_left, concentration])
# In[718]:
ls = ScipyDirect({})
nls_status = IndexedStruct()
nls = Newton({'is_linear' : True}, lin_solver=ls, status=nls_status)
tss = SimpleTimeSteppingSolver({'t0': 0.0, 't1': 50.0, 'n_step': 50}, nls=nls, context=pb, verbose=True, step_hook=pb.update_time_stepper)
pb.set_bcs(ebcs=ebcs)
pb.set_solver(tss)
pb.time_update(tss.ts)
pb.set_ics(Conditions([ic]))
pb.update_materials()
# In[719]:
pb.set_bcs(ebcs=ebcs)
pb.set_solver(nls)
pb.set_solver(tss)
pb.time_update(ebcs=Conditions(ebcs=ebcs))
status = IndexedStruct()
variables = pb.solve(status=status)
print('Nonlinear solver status:\n', status) # Assuming nls_status is not defined, using 'status' instead
print('Stationary solver status:\n', status) # Assuming nls_status is not defined, using 'status' instead
pb.save_state('Concentration.vtk', variables)
Thanking you
Himalaya singh
