plot Expression not correct

Question

I use dolfin 1.2.0 (Python-Version) and Ubuntu 12.04. I'm confused with the plot(Expression,mesh). I want to plot a Expression to see if my function or values is right or not, but got always the wrong values graphic. I paste my short code here. If someone can help, I'll be very appreciated!!!

from dolfin import *
import numpy

#Create mesh
mesh = Mesh("model.xml")
subdomains = MeshFunction("size_t", mesh, "model_physical_region.xml")

#Define function space with Nedelec (edge) elements
V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 2)

#Define test and trial functions
v = TestFunction(X_h)
u = TrialFunction(X_h)

# Define the current values in the coils
class Current_Coils_SubDomain(Expression):  
     def eval_cell(self, values, x, ufc_cell):
          j = 1.0
          i = 0     
          region_index = subdomains[ufc_cell.index]
          x,y,z = [x[0], x[1], x[2]]                    
          if region_index == 1 or region_index == 2:                                
             r = pow(x**2 + y**2, 0.5)
             values[0] = -y/r*j
             values[1] = x/r*j
             values[2] = 0.0
          else:
             values[0] = 0.0
             values[1] = 0.0
             values[2] = 0.0                    
     def value_shape(self):
          return(3,)

#Instantiate the current
J = Current_Coils_SubDomain()

# Test the current
#print J, J.str(), J.ufl_element(), J.value_shape()
plot(J, mesh, interactive=True)

And what I got is the picture like this one:

And I don't know how to correct it, my xml data was drawed and meshed with Gmsh(2.5.1) and converted with "dolfin-convert". The values I want to set are in the round coils and they should be around over the x-y achse.

Comments

Hi, could you please attach the plot and your xml data?

---

Hi Mirok,
I don't know how to upload attachment here. I can just put my xml-data here.
1. model.xml
enter link description here

2.model_physical_region.xml
enter link description here

3.the plot file(sorry, I changed some mesh, so now it gets worse....It shows now all with value 0.)
enter link description here

---

I get "access denied" when following either of the links.

---

Do you have email? I can send you per mail.
Because I save the files in the bitbucket, those are link of files. I think people must have a account of that website "bitbucket.org".

---

Hi Mirok,
I think i find you in the bitbucket and add you to my user list. Can you try one more?

Answer 1

I may be wrong, but I believe that adding the code

f = Function(V)
f.interpolate(J)
plot(f)

may help. The Expression object needs to be projected onto a function space (this may be incorrect terminology).

Let me know if this helps.

Comments

I have tried and by my side it is still failed....Anyway thanks for your suggestion.

Answer 2

Hi, I think the problem with your code is that during computation of the expansion coefficients the same dof is visited as part of the ufc_cell with index 1 or 2 but then it is also visited as part of cell with index 3, 4 and therefore the values are set to zeros. This is quite plausible if your mesh for the coils is such that every tetrahedron in the coils has at least one member in the region 3 or 4. I checked your mesh and that seems to be the case. So the first suggestion is to refine the coil mesh such that there tets only surrounded by other tets with indices 1, 2. The other suggestion which works with the current mesh is:

from dolfin import *
import numpy

#Create mesh
mesh = Mesh("model.xml")
subdomains = MeshFunction("size_t", mesh, "model_physical_region.xml")

# Define the current values in the coils
class Current_Coils_SubDomain(Expression):  
     def eval_cell(self, values, x, ufc_cell):
          j = 1.0
          i = 0     
          region_index = subdomains[ufc_cell.index]
          x,y,z = [x[0], x[1], x[2]]                    
          if region_index == 1 or region_index == 2:                                
             r = pow(x**2 + y**2, 0.5)
             values[0] = -y/r*j
             values[1] = x/r*j
             values[2] = 0.0
          else:
             values[0] = 0.0
             values[1] = 0.0
             values[2] = 0.0                    
     def value_shape(self):
          return(3,)

#Instantiate the current
J = Current_Coils_SubDomain()

V = VectorFunctionSpace(mesh, "DG", 0)
u = interpolate(J, V)
W = VectorFunctionSpace(mesh, "CG", 1)
v = project(u, W)
plot(v, interactive=True) 

Comments

Hi Mirok,
I'll try your both suggestion and later tell you the result. Thanks for your help again.

---

I think that at moment I can't find a better way except yours to solve this question. Can I ask a more question about "out of memory"?
I used the same xml files to solve a linear variational system and get the error like:

UMFPACK V5.4.0 (May 20, 2009): ERROR: out of memory
[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Error in external library!
[0]PETSC ERROR: umfpack_UMF_numeric failed!
[0]PETSC ERROR: ------------------------------------------------------------------------

Have you met the same problem beore?

---

Hi, your system might be too big for direct solver. See these threads 1, 2.