esys-escript

problem with setValue

Asked by Louise Olsen

I would like to use setTaggedValue to define the boundary conditions in my PDE. However I get an error message if I use these commands:

mypde = LinearPDE(mydomain)

# setting boundary conditions

boundary_cond=Vector(0.,Function(mydomain))
boundary_cond.setTaggedValue("Top",[0.,0.,1.])
boundary_cond.setTaggedValue("Bottom",[0.,0.,1.])
msk_point=whereZero(mesh[2])*whereZero(mesh[1])*whereZero(mesh[0])*[1.,1.,1.]
boundary_cond+=msk_point

mypde.setValue(q=boundary_cond)

This is the error message:

 run-escript NoDamageCylinderNoFlaw.py
Traceback (most recent call last):
  File "NoDamageCylinderNoFlaw.py", line 84, in <module>
    mypde.setValue(q=boundary_cond)
  File "escriptcore/py_src/linearPDEs.py", line 2053, in setValue
  File "escriptcore/py_src/linearPDEs.py", line 1298, in setValue
esys.escriptcore.linearPDEs.IllegalCoefficientFunctionSpace: Coefficient q:Unable to interpolate coefficient to function space Finley_DegreesOfFreedom on FinleyMesh

However the code below works:

msk_bottom=whereZero(mesh[2])*[0.,0.,1.]
msk_top=whereZero(mesh[2]-lz)*[0.,0.,1.]
msk_point=whereZero(mesh[2])*whereZero(mesh[1])*whereZero(mesh[0])*[1.,1.,1.]
boundary_cond = msk_bottom + msk_top + msk_point

Thanks, Louise

Question information

Language:
English Edit question
Status:
Solved
For:
esys-escript Edit question
Assignee:
No assignee Edit question
Solved by:
Bob
Solved:
Last query:
Last reply:
Revision history for this message
Best Bob (caltinay) said : 		#1

The coefficient q needs to be defined on the domain nodes (degrees of freedom) rather than the elements (=Function).
To make your first example work, you only need to change the definition of boundary_cond
from:
boundary_cond=Vector(0.,Function(mydomain))
to:
boundary_cond=Vector(0., ContinuousFunction(mydomain))

In your second example boundary_cond is implicitly defined on ContinuousFunction because msk_bottom, msk_top and msk_point are defined via mesh[i] which I presume is domain.getX().

Revision history for this message
Louise Olsen (l-kettle1) said : 		#2

Thanks Cihan Altinay, that solved my question.

Revision history for this message
Joel Fenwick (j-fenwick1) said : 		#3

In the next release, if you print out the FunctionSpace of a data object, it will give an indication of how the relevant FunctionSpace was made.

eg:
  print (d.getFunctionSpace())
  Ripley_Nodes [ContinuousFunction(domain)] on ripley::Rectangle

(as I say: next release).