Defining spatially variable scalars using Expressions for 2D and 3D problems

Question

Hello All,

I am trying to define a material property that varies in space (for my problem I have an inclusion and a matrix, each with its own material property on a unit cube mesh). I tried using the Expressionclass as follows:

class materialPoisson(Expression):

  def __init__(self, a_NuMatrix, a_NuInclusion):

    self.m_NuMatrix     = a_NuMatrix
    self.m_NuInclusion  = a_NuInclusion

  def eval(self, value, x):

    if myInclusion.inside(x, on_boundary):
      Nu0 = self.m_NuInclusion
    else:
      Nu0 = self.m_NuMatrix

    value = Nu0

  def value_shape(self):
    return (1,)

The function myInclusion.inside() is just checking whether the point lies inside the inclusion domain or outside of it. I want to now use this in defining the variational form, and try to interpolate it as follows:

Nu = interpolate(Nu, V)

where

 V = VectorFunctionSpace(mesh, "CG", 1)

I end up always getting the following error:

*** -------------------------------------------------------------------------
*** Error:   Unable to interpolate function into function space.
*** Reason:  Dimension 0 of function (1) does not match dimension 0 of function space (3).
*** Where:   This error was encountered inside FunctionSpace.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  3b6582dfb45139c906c13b9ad57395632a2090f4
*** -------------------------------------------------------------------------

How do I interpolate the scalar expression on a vector function space then? Or should I be using a different approach to defining a spatially varying scalar field which I can use for defining material properties (say) in variational forms?

Any tips or pointers will be much appreciated. :)

Thanks

Comments

Hi, if it okay to consider your material property as piecewise constant then this answer might be helpful

---

If Nu is a scalar (according to the shape you have defined), you cannot interpolate it to a VectorFunctionSpace.

---

Thanks MiroK - I followed the procedure on that link and I have the following code now:

class materialLambda(Expression):

def __init__(self, a_Inclusion, a_LambdaMatrix, a_LambdaInclusion):

  self.m_Inclusion        = a_Inclusion
  self.m_LambdaMatrix     = a_LambdaMatrix
  self.m_LambdaInclusion  = a_LambdaInclusion

def eval(self, value, x):

  if isinstance(self.m_Inclusion, Inclusion):

    if self.m_Inclusion.inside(x):
      Lambda = self.m_LambdaInclusion
    else:
      Lambda = self.m_LambdaMatrix

  elif isinstance(self.m_Inclusion, Superquadric):

    if self.m_Inclusion.insideFunction(x[0],x[1],x[2]):
      Lambda = self.m_LambdaInclusion
    else:
      Lambda = self.m_LambdaMatrix

  value = Lambda

where the classes Inclusion and Superquadric are SubDomains that have a explicitly defined inside function. I now interpolate this as:

D = FunctionSpace(mesh,"DG",0)
Lambdaf = materialLambda(myInclusion, LambdaMatix, LambdaInclusion)
Lambda = interpolate(Lambdaf, D)

However I think this interpolation just gives me 0 everywhere. I plot this Lambda using plot(Lambda) and get zeros everywhere. Is there something else I am missing ?

---

Hi, in eval you need to have value[0]=Lambda, see here for more. Sorry that I didn't spot this in your first code.

---

Hello MiroK

Thanks !! That solved the issue I was having - I can now set this up to interpolate a piecewise constant on my domain.