jorgensd / dolfinx_mpc

Extension for dolfinx to handle multi-point constraints.
https://jorgensd.github.io/dolfinx_mpc/
MIT License
31 stars 12 forks source link

Possibility of periodic boundary condition for cylindrical polar coordinates #19

Open tbrandvik opened 2 years ago

tbrandvik commented 2 years ago

I am trying to implement a periodic boundary condition for the 3D elastic equation. The only difference compared to the existing example in the demos is that the boundaries are periodic in cylindrical polar coordinates (e.g. the slave nodes are at the same coordinates as the master nodes, but rotated through an angle). This leads to the following kind of relationship between the displacements of the periodic nodes, assuming the second node is the first node rotated about the x-axis:

$$ \begin{align} dy_2 &= a\cdot dy_1 + b\cdot dz_1 \ dz_2 &= c\cdot dy_1 + d\cdot dz_1 \end{align} $$

The coefficients (a, b, c, d) would be fixed and set by the rotation matrix.

I was wondering if anyone has any ideas on how to implement this with dolfinx_mpc? Any help would be greatly appreciated!

jorgensd commented 2 years ago

You should be able to do this with https://github.com/jorgensd/dolfinx_mpc/blob/master/python/dolfinx_mpc/multipointconstraint.py#L97 where you use the relation function to map (x1,y1) to radial scaling Rotation_matrix (x1,y1).

bleyerj commented 2 years ago

I think there might be a confusion. Indeed one can map (x1, y1) to (x2, y2) through some rotation but one should also enforce

(ux(x2, y2), uy(x2,y2)) = Rotation_matrix*(ux(x1, y1), uy(x1, y1))

I don't see how this could be achieved with the current functionalities since the dof corresponding to one subspace should be a linear combination of the subspace dofs at the opposite point

jorgensd commented 2 years ago

That is not a periodic condition then, and you would need to make a custom function to do this (using for instance https://github.com/jorgensd/dolfinx_mpc/blob/master/python/dolfinx_mpc/multipointconstraint.py#L217-L245). A periodic condition (in my opinion) is:

ux (x1, y1) = ux(Rotation_matrix * (x1,y1))
uy (x1, y1) = uy(Rotation_matrix * (x1,y1))
bleyerj commented 2 years ago

Alright, it is periodic with respect to the local tangent frame components, not the global ones.

But so with https://github.com/jorgensd/dolfinx_mpc/blob/master/python/dolfinx_mpc/multipointconstraint.py#L217-L245), I don't see a way of enforcing for one subspace_slave (e.g V.sub(0)) at point numpy.array([x2, y2], dtype=numpy.float64).tobytes() a linear combination of two subspace_master (e.g. V.sub(0) and V.sub(1)) at the same point numpy.array([x1, y1], dtype=numpy.float64).tobytes() Am I missing something ?

jorgensd commented 2 years ago

The subspace_slave and subspace_master is just for convenience for cases where you can isolate the components. If you do not set these, it will work with dofs on the full space (dofs unrolled with block size).

bleyerj commented 2 years ago

Got it. But then the coefficient alpha of the MPC equation is only a scalar or it can be a matrix with dimension block_size x block_size ?

jorgensd commented 2 years ago

You would need to code each coefficient for each degree of freedom separately, i.e. for each (x1, y1) we have a slave dof_0(x1, y1)_x, masters dof_1(map(x1, y1))_x, dof_1(map(x1,y1))_y, we would have coefficients rotation_matrix[0,0], rotation_matrix[1, 0] (if I've not missed something).

bleyerj commented 2 years ago

So I would set two MPCs for the same slave dof_0(x1,y1)_x :

would it "sum" both master contributions ?

jorgensd commented 2 years ago

Maybe the create_general_constraint is not the right abstraction for your problem (I need to think about it) and revisit the code.

The most general way of adding a constraint is by using flat arrays with offsets indicating what dofs relate to another set of dofs: https://github.com/jorgensd/dolfinx_mpc/blob/master/python/dolfinx_mpc/multipointconstraint.py#L39-L61 which is easy in serial, but requires your own MPI communication to determine relationships in parallel.

tbrandvik commented 2 years ago

Thank you both for looking into this! I think I can see how it might be done using flat arrays so I will have a go at that for now.