Hi Phil,
On 04/07/2014 09:33 AM, seismo_phil wrote:
Hello,
I'm still working on the same thermal diffusion problem I referred to in the "thermal diffusion - Dirichlet vs. Neumann BCs" thread. The problem is attached file "flndrs2d.py" and *appears* to give a reasonable solution, but I still need to calculate heat flux out the top of the model. So, following Robert's suggestion i added a post_process hook that evaluates d_surface_flux:
def post_process(out, problem, state, extend=False):
"""
Calculate gradient of the solution.
"""
flux = problem.evaluate('d_surface_flux.2.Gamma_Top(coef.val, t)',
mode='eval', copy_materials=False, verbose=False)
print flux
When I call this using my rather complicated materials definition in the attached file, d_surface flux fails with the following error:
File "/Users/philcummins/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/sfepy/discrete/functions.py", line 114, in get_constants
nqp = qps.shape[ig][1]KeyError: 1
(after a lengthy traceback). If I change the materials definition to a very
There is a bug in ConstantFunctionByRegion that occurs when evaluating a term in a region that is not listed in the material definition, although the region is a subset of the regions in the list (so that it should work).
simple one:
materials = {
'flux' : ({'val' : 0.03},), 'coef' : ({'val' : 1.0},),}
Then d_surface_flux appears to evaluate OK, but it just returns a single value of zero. (I am expecting a vector). I'm not sure what to try next, any suggestions?
The easiest way (before the bug is fixed) is to use a function for defining the coefficients.
def get_coef(ts, coors, mode=None, **kwargs): if mode != 'qp': return
indx = kwargs['group_indx']
term = kwargs['term']
regions = kwargs['problem'].domain.regions
val = np.empty((coors.shape[0], 1, 1), dtype=np.float64)
print term.region.name
for key, par in prms.iteritems():
region = regions[key]
for ig in term.region.igs:
if ig not in region.igs: continue
print region.name, indx[ig]
val[indx[ig]] = par['k'] / (par['rho'] * par['cp'])
return {'val' : val, 'valI' : val * np.eye(2)}I am attaching my modifications of your script that seems to work. Note also:
- pass mode='el' to problem.evaluate() to get flux components for each face
- d_surface_flux requires its material parameter to be a tensor, so I defined 'valI' in the function above.
Let us know if that works...
Cheers, r.
Many thanks
- Phil
P.S. I have also attached a coarsened mesh file, in case anyone wants to try to run it. But you need to remove the "2km" in the name of the mesh file in flndrs2d.py if you want it to work.