This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Regarding structure of VectorFunctionSpace

+3 votes

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?

asked Jul 21, 2014 by wensby FEniCS Novice (350 points)

3 Answers

+3 votes

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.

answered Jul 23, 2014 by Garth N. Wells FEniCS Expert (35,930 points)

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.

0 votes

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.

answered Jul 25, 2014 by wensby FEniCS Novice (350 points)

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 :)

+1 vote

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.

answered Jul 25, 2014 by Charles FEniCS User (4,220 points)
...