Using the dw_diffusion_coupling term with inhomogenous materials

Hello,
I've been using SfePy for about a week - it's a great piece of work.
I'm working on a 2D problem involving a dw_diffusion_coupling termhttp://sfepy.org/doc-devel/src/sfepy/terms/termsLaplace.html#sfepy.terms.termsLaplace.DiffusionCoupling. For a homogenous material everything works as expected using the material definition
'mat_diffusion_coupling': ({'f': np.array([[0], [0]])},),
For an inhomogenous material, since the shape of this local parameter is (2, 1), I defined a function which takes the array of coordinates - coors - and returns an array of shape (coors.shape[0], 2, 1) like so
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs):
return {'f': np.zeros((coors.shape[0], 2, 1))}
This gives the error
File ".../sfepy/terms/termsLaplace.py", line 199, in dw_fun
status = terms.mulATB_integrate(out, vg.bfg, mat * val, vg)
ValueError: operands could not be broadcast together with shapes (5151,2,1) (5000,8,1,1)
, so I expect I'm doing something silly! How do I get this to work?
(I've attached a minimal version of the code giving the error.)
Many thanks, Ben

Hi Ben,
for an inhomogenous material the shape of the material parameter should be (dim, dim), it corresponds to the permeability tensor K_ij. So, in your code
'mat_diffusion_coupling': ({'f': np.eye(2)},),
or
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs): return {'f': np.zeros((coors.shape[0], 2, 2))}
will work.
Regards, Vladimir
On 03/12/2013 04:40 PM, Ben Derrett wrote:
Hello,
I've been using SfePy for about a week - it's a great piece of work.
I'm working on a 2D problem involving a dw_diffusion_coupling term http://sfepy.org/doc-devel/src/sfepy/terms/termsLaplace.html#sfepy.terms.termsLaplace.DiffusionCoupling. For a homogenous material everything works as expected using the material definition
'mat_diffusion_coupling': ({'f': np.array([[0], [0]])},),
For an inhomogenous material, since the shape of this local parameter is (2, 1), I defined a function which takes the array of coordinates - coors - and returns an array of shape (coors.shape[0], 2, 1) like so
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs): return {'f': np.zeros((coors.shape[0], 2, 1))}
This gives the error
File ".../sfepy/terms/termsLaplace.py", line 199, in dw_fun status = terms.mulATB_integrate(out, vg.bfg, mat * val, vg) ValueError: operands could not be broadcast together with shapes (5151,2,1) (5000,8,1,1)
, so I expect I'm doing something silly! How do I get this to work?
(I've attached a minimal version of the code giving the error.)
Many thanks, Ben
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sfepy-devel...@googlegroups.com. To post to this group, send email to sfepy...@googlegroups.com. Visit this group at http://groups.google.com/group/sfepy-devel?hl=en. For more options, visit https://groups.google.com/groups/opt_out.

Hi Vladimir,
Thank you for replying. According to the documentation the dw_diffusion_coupling term represents the integral
\int_\Omega pK_j\nabla_j q,
with the summation convention. So it states that K is a vector, rather than a matrix. If this is wrong, which integral does the dw_diffusion_coupling term represent?
Thanks again, Ben
On Tuesday, March 12, 2013 5:24:41 PM UTC, vladimir lukes wrote:
Hi Ben,
for an inhomogenous material the shape of the material parameter should be (dim, dim), it corresponds to the permeability tensor K_ij. So, in your code
'mat_diffusion_coupling': ({'f': np.eye(2)},),
or
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs): return {'f': np.zeros((coors.shape[0], 2, 2))}
will work.
Regards, Vladimir
On 03/12/2013 04:40 PM, Ben Derrett wrote:
Hello,
I've been using SfePy for about a week - it's a great piece of work.
I'm working on a 2D problem involving a dw_diffusion_coupling term <
http://sfepy.org/doc-devel/src/sfepy/terms/termsLaplace.html#sfepy.terms.ter....
For a homogenous material everything works as expected using the material definition
'mat_diffusion_coupling': ({'f': np.array([[0], [0]])},),
For an inhomogenous material, since the shape of this local parameter is (2, 1), I defined a function which takes the array of coordinates - coors - and returns an array of shape (coors.shape[0], 2, 1) like so
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs): return {'f': np.zeros((coors.shape[0], 2, 1))}
This gives the error
File ".../sfepy/terms/termsLaplace.py", line 199, in dw_fun status = terms.mulATB_integrate(out, vg.bfg, mat * val, vg) ValueError: operands could not be broadcast together with shapes (5151,2,1) (5000,8,1,1)
, so I expect I'm doing something silly! How do I get this to work?
(I've attached a minimal version of the code giving the error.)
Many thanks, Ben
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sfepy-d...@googlegroups.com javascript:. To post to this group, send email to sfep...@googlegroups.comjavascript:.
Visit this group at http://groups.google.com/group/sfepy-devel?hl=en. For more options, visit https://groups.google.com/groups/opt_out.

Hi Ben,
I apologize, I was wrong, in dw_diffusion_coupling term a vector, not a matrix, must be used as the material parameter.
I see the problem in the material function, try this:
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs): if mode == 'qp': return {'f': np.zeros((coors.shape[0], 2, 1))}
else:
return None
The values are returned only in 'quadrature points' mode.
Regards Vladimir
On 03/12/2013 08:21 PM, Ben Derrett wrote:
Hi Vladimir,
Thank you for replying. According to the documentation the dw_diffusion_coupling term represents the integral
\int_\Omega pK_j\nabla_j q,
with the summation convention. So it states that K is a vector, rather than a matrix. If this is wrong, which integral does the dw_diffusion_coupling term represent?
Thanks again, Ben
On Tuesday, March 12, 2013 5:24:41 PM UTC, vladimir lukes wrote:
Hi Ben, for an inhomogenous material the shape of the material parameter should be (dim, dim), it corresponds to the permeability tensor K_ij. So, in your code 'mat_diffusion_coupling': ({'f': np.eye(2)},), or def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs): return {'f': np.zeros((coors.shape[0], 2, 2))} will work. Regards, Vladimir On 03/12/2013 04:40 PM, Ben Derrett wrote: > Hello, > > I've been using SfePy for about a week - it's a great piece of work. > > I'm working on a 2D problem involving a dw_diffusion_coupling term > <http://sfepy.org/doc-devel/src/sfepy/terms/termsLaplace.html#sfepy.terms.termsLaplace.DiffusionCoupling <http://sfepy.org/doc-devel/src/sfepy/terms/termsLaplace.html#sfepy.terms.termsLaplace.DiffusionCoupling>>. > For a homogenous material everything works as expected using the > material definition > > 'mat_diffusion_coupling': ({'f': np.array([[0], [0]])},), > > > For an inhomogenous material, since the shape of this local parameter is > (2, 1), I defined a function which takes the array of coordinates - > coors - and returns an array of shape (coors.shape[0], 2, 1) like so > > def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs): > return {'f': np.zeros((coors.shape[0], 2, 1))} > > > This gives the error > > File ".../sfepy/terms/termsLaplace.py", line 199, in dw_fun > status = terms.mulATB_integrate(out, vg.bfg, mat * val, vg) > ValueError: operands could not be broadcast together with shapes > (5151,2,1) (5000,8,1,1) > > > , so I expect I'm doing something silly! How do I get this to work? > > (I've attached a minimal version of the code giving the error.) > > Many thanks, > Ben > > -- > You received this message because you are subscribed to the Google > Groups "sfepy-devel" group. > To unsubscribe from this group and stop receiving emails from it, send > an email to sfepy-d...@googlegroups.com <javascript:>. > To post to this group, send email to sfep...@googlegroups.com <javascript:>. > Visit this group at http://groups.google.com/group/sfepy-devel?hl=en <http://groups.google.com/group/sfepy-devel?hl=en>. > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>. > >
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sfepy-devel...@googlegroups.com. To post to this group, send email to sfepy...@googlegroups.com. Visit this group at http://groups.google.com/group/sfepy-devel?hl=en. For more options, visit https://groups.google.com/groups/opt_out.

Hi Vladimir,
On 03/12/2013 09:47 PM, Vladimír Lukeš wrote:
Hi Ben,
I apologize, I was wrong, in dw_diffusion_coupling term a vector, not a matrix, must be used as the material parameter.
I see the problem in the material function, try this:
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs): if mode == 'qp': return {'f': np.zeros((coors.shape[0], 2, 1))}
else: return None
The values are returned only in 'quadrature points' mode.
Just a small note - the "return None" branch does not need to be there explicitly. Usually, if only 'qp' mode is used in a material function, I use the pattern:
if mode != 'qp': return
...
Also, Ben, the call modes are described in [1], as well as the syntax of other function hooks.
Cheers, r.

On 03/12/2013 04:40 PM, Ben Derrett wrote:
Hello,
I've been using SfePy for about a week - it's a great piece of work.
I'm working on a 2D problem involving a dw_diffusion_coupling termhttp://sfepy.org/doc-devel/src/sfepy/terms/termsLaplace.html#sfepy.terms.termsLaplace.DiffusionCoupling. For a homogenous material everything works as expected using the material definition
'mat_diffusion_coupling': ({'f': np.array([[0], [0]])},),
For an inhomogenous material, since the shape of this local parameter is (2, 1), I defined a function which takes the array of coordinates - coors - and returns an array of shape (coors.shape[0], 2, 1) like so
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs):
return {'f': np.zeros((coors.shape[0], 2, 1))}
Another thing: the attached code is just a minimal example without realistic BCs etc., right? If I run the attached file (with the mode correction posted by Vladimir), it does nothing because the loads as well as material parameters are zero.
r.

Thank you both! It works well now.
Yes, Robert. I cut out every extraneous feature of my problem.
Ben
On Tuesday, March 12, 2013 9:20:48 PM UTC, Robert Cimrman wrote:
On 03/12/2013 04:40 PM, Ben Derrett wrote:
Hello,
I've been using SfePy for about a week - it's a great piece of work.
I'm working on a 2D problem involving a dw_diffusion_coupling term<
http://sfepy.org/doc-devel/src/sfepy/terms/termsLaplace.html#sfepy.terms.ter....
For a homogenous material everything works as expected using the
material
definition
'mat_diffusion_coupling': ({'f': np.array([[0], [0]])},),
For an inhomogenous material, since the shape of this local parameter is (2, 1), I defined a function which takes the array of coordinates -
coors -
and returns an array of shape (coors.shape[0], 2, 1) like so
def diffusion_coupling_coefficient(ts, coors, mode=None, **kwargs):
return {'f': np.zeros((coors.shape[0], 2, 1))}
Another thing: the attached code is just a minimal example without realistic BCs etc., right? If I run the attached file (with the mode correction posted by Vladimir), it does nothing because the loads as well as material parameters are zero.
r.
participants (3)
-
Ben Derrett
-
Robert Cimrman
-
Vladimír Lukeš