r""" An adaptation to the Transient Laplace equation. The same example as time_poisson.py, but using the short syntax of keywords, and explicit time-stepping. Find :math:T(t) for :math:t \in [0, t_{\rm final}] such that: .. math:: \int_{\Omega} s \pdiff{c}{t} + \int_{\Omega} D \nabla s \cdot \nabla c = \int_{\Gamma_N} s g \;, \quad \forall s \;. where :math:g is the given flux, :math:g = \nabla T \cdot \ul{n}. See the tutorial section :ref:poisson-weak-form-tutorial for a detailed explanation. View the results using:: \$ ./postproc.py -b -d't,plot_warp_scalar,rel_scaling=1' --wireframe --view=-90,90,1.5,0,0,0 --roll=0 laplace_1d.h5 """ from __future__ import absolute_import from sfepy import data_dir import numpy as nm from sfepy.discrete.fem import Mesh from sfepy.discrete.fem.meshio import UserMeshIO def mesh_hook(mesh, mode): """ Generate the 1D mesh. """ if mode == 'read': n_nod = 201 #n_nod = 101 # n_nod = 51 coors = nm.linspace(0.0, 200.0, n_nod).reshape((n_nod, 1)) #coors = nm.linspace(0.0, 100.0, n_nod).reshape((n_nod, 1)) # coors = nm.linspace(0.0, 50.0, n_nod).reshape((n_nod, 1)) conn = nm.arange(n_nod, dtype=nm.int32).repeat(2)[1:-1].reshape((-1, 2)) mat_ids = nm.zeros(n_nod - 1, dtype=nm.int32) descs = ['1_2'] mesh = Mesh.from_data('laplace_1d', coors, None, [conn], [mat_ids], descs) return mesh elif mode == 'write': pass def get_ic(coor, ic): """Non-constant initial condition.""" import numpy as nm ##arbitrary scale factor scalefactor = 10 ##define gaussian magnitude scaled to cmax a = 1.0 C = 30*a A = 1/(C*nm.sqrt(2*nm.pi)) B = 0 cinit = nm.zeros(len(coor)) for i in range(len(coor)): cinit[i] = scalefactor*A*nm.exp(-(coor[i]-B)*(coor[i]-B)/(2*C*C)) # print (cinit) return cinit #filename_mesh = data_dir + '/meshes/3d/cylinder.mesh' filename_mesh = UserMeshIO(mesh_hook) materials = { 'coef' : ({'D' : 1.0},), 'LoadRHS' : ({'.val' : 0.0},), 'LoadLHS' : ({'.val' : 0.0},), } regions = { 'Omega' : 'all', 'Gamma_Left' : ('vertices in (x < 0.00001)', 'vertex'), 'Gamma_Right' : ('vertices in (x > 199.99999)', 'vertex'), # 'Gamma_Right' : ('vertices in (x > 49.99999)', 'vertex'), # 'Gamma_Left' : ('vertex '), # 'Gamma_Right' : ('vertex 101'), } fields = { 'concentration' : ('real', 1, 'Omega', 1), } variables = { 'c' : ('unknown field', 'concentration', 0, 1), 's' : ('test field', 'concentration', 'c'), } # ebcs = { # 'c1' : ('Gamma_Left', {'c.0' : 0.0}), # 'c2' : ('Gamma_Right', {'c.0' : 0.0}), # } ics = { 'ic' : ('Omega', {'c.0' : 'get_ic'}), } functions = { 'get_ic' : (get_ic,), } integrals = { 'i' : 1, } equations = { 'concentration' : """dw_volume_dot.i.Omega( s, dc/dt ) + dw_laplace.i.Omega( coef.D, s, c ) = dw_point_load.0.Gamma_Right( LoadRHS.val, s) + dw_point_load.0.Gamma_Left( LoadLHS.val, s) """ ### for 1D 0 flux BCs # dw_point_load.0.Gamma_Right( Load.val, s) + dw_point_load.0.Gamma_Left( Load.val, s) ### for 1D # = dw_point_load.0.Gamma_Right( Load.val, s) ### for 2D/3D Meshes # = dw_surface_integrate.2.Gamma_Right( flux.val, s ) } solvers = { 'ls' : ('ls.scipy_direct', {}), 'newton' : ('nls.newton', { 'i_max' : 1, 'is_linear' : True, }), 'ts' : ('ts.simple', { 't0' : 0.0, 't1' : 300.0, 'dt' : 5.0, 'n_step' : None, }), } output_dir = './output' # set this to a valid directory you have write access to options = { 'ls' : 'ls', 'ts' : 'ts', 'save_steps' : 2, 'output_format' : 'vtk', # string, output directory 'output_dir' : output_dir, # 'step_hook' : '', }