Regarding structure of VectorFunctionSpace

Question

I have two simple subdomains representing two separate solid mechanical bodies and I'm trying to set different body force parameters (based on gravity and density) for each subdomain. Right now I'm looping through each cell and setting the x-, y- and z-component of the body force for each cell based on the subdomains density parameter which I've set earlier.

# Subdomain Body Forces
VFS = VectorFunctionSpace(mesh, 'DG', 0)
B = Function(VFS)
cell_count = len(subdomains.array())
for cell_no in range(cell_count):
   # Fetch subdomain number
   subdomain_no = subdomains.array()[cell_no]
   # Set x-, y- and z-component of the cell
   B.vector()[cell_no + (0 * cell_count)] = 0 * rho.vector()[cell_no]
   B.vector()[cell_no + (1 * cell_count)] = 0 * rho.vector()[cell_no] 
   B.vector()[cell_no + (2 * cell_count)] = -9.8 * rho.vector()[cell_no]

I was told that the B.vector() would be structured as:

[x1 x2 ... xn y1 y2 ... yn z1 z2 ... zn]

But perhaps this is not the case. When I do like this, only one of my bodies receives a force, strangely pointing down and to the back as if the body forces of the subdomains turned out

B(subdomain0) = (0, 0, 0) (times the cell-specific density)
B(subdomain1) = (0, -9.8, -9.8) (times the cell-specific density)

Any suggestions on where I'm thinking incorrectly? How is B.vector() structured?

Answer 1

You cannot make any assumptions on how the degrees of freedom are ordered. You need to use the DofMap if you want to extract particular values.

Comments

Thank you! So you are saying that if I just want to assign a scalar using a Function (as they do here), it is safe to assume the structure of the function's vector, but with VectorFunction, all bets are off?

All I really want to do is have different B = Constant((x, y, z)) for the different subdomains. I don't want to extract particular values, but rather set particular values B = (x, y, z) for each subdomain (I guess subsequently for every cell).

But if I still need to use DofMap, I have to admit I'm not sure on how to. I seem to find very little examples on it. Could this approach work in my situation too?

---

No, I said 'You cannot make any assumptions on how the degrees of freedom are ordered'. This question has been asked numerous times before on this forum. Please search for past questions.

Answer 2

Although one shouldn't assume the structure of B.vector(), I nevertheless solved this problem when assigning values as if the x-, y- and z-components for a specific cell were right next to each other:

for cell_no in range(cell_count):
   # Fetch subdomain number
   subdomain_no = subdomains.array()[cell_no]
   # Set x-, y- and z-component of the cell
   B.vector()[(3 * cell_no) + 0] =  0.0 * rho.vector()[cell_no]
   B.vector()[(3 * cell_no) + 1] =  0.0 * rho.vector()[cell_no]
   B.vector()[(3 * cell_no) + 2] = -9.8 * rho.vector()[cell_no]

Of course, this may very likely be a solution unique to my particular program and not work in the general case.

Comments

This is a very volatile solution...

---

I see that now, so I would not recommend this solution. I wouldn't mind it being down-voted.

To be honest, this was only initially a problem since I wanted to see if I could have different gravity vectors for each subdomain. But since then I settled with a uniform gravity vector over all subdomains. Your answer works exemplary in that case!

---

You can handle a non-uniform gravity field by creating an expression for g, interpolating that onto a function and doing B = rho * g from there as well :)

Answer 3

I'm guessing by 9.8 times rho you are getting a body force...

Create a vector constant for body acceleration and multiply it by rho, then assign that to B (I'm guess for body force)

// use a vector to automatically handle your concerns with xyz
Constant g(0,0,-9.8);

// use the operators for the Functions to perform the vector operation
B = rho * g;

Now your code isn't susceptible to internal changes in the project.