Hi!
Actually, I would like to compute for each node of the mesh the spherical coordinates and associated local spherical basis. As I am in dynamic case, I would like them to be updated at each solver iteration.
Here are my questions:
- I understood, but maybe I am wrong, that I needed to define the nodal spherical coordinates and basis as finite element functions to be sure they would be updated. Is it correct? Is there another way?
- In the below code, I cannot associate a ‘fe.as_vector’ (dimension 3) to ‘e_r_vect’ (np.array([0.,0.,0…]) ) and get a ‘setting an array element with a sequence’ ValueError. As e_r_vect is supposed to contain the 3d components for all node (I thought shape would be : (number of node, 3) ), how is it possible to access them precisely from the shape(3 * number of nodes,)-structure of ‘e_r_vect’?
For now, I tried the following code:
Define vector space and associated spherical basis vectors functions
VectorSpace_DG1 = VectorFunctionSpace(mesh, “Lagrange”, 1)
e_r = Function(VectorSpace_DG1)
e_theta = Function(VectorSpace_DG1)
e_phi = Function(VectorSpace_DG1)
Define scalar space and associated functions
ScalarSpace_DG1 = FunctionSpace(mesh, “Lagrange”, 1)
r = Function(ScalarSpace_DG1)
theta = Function(ScalarSpace_DG1)
phi = Function(ScalarSpace_DG1)
coordinates_all_nodes = np.array([ np.array([node.point().x(), node.point().y(), node.point().z()]) for dof, node in enumerate(fe.vertices(mesh))])
Dynamic spherical coordinates?
r_vec = r.vector()
theta_vec = theta.vector()
phi_vec = phi.vector()
for dof, node in enumerate(vertices(fenics_mesh)):
r_vec[dof] = sqrt( coordinates_all_nodes[dof][0]*coordinates_all_nodes[dof][0] +
coordinates_all_nodes[dof][1]*coordinates_all_nodes[dof][1] +
coordinates_all_nodes[dof][2]*coordinates_all_nodes[dof][2] )
theta_vec[dof] = acos(coordinates_all_nodes[dof][2]/r_vec[dof])
phi_vec[dof] = atan(coordinates_all_nodes[dof][1]/coordinates_all_nodes[dof][0])
r.vector()[:] = r_vec
theta.vector()[:] = theta_vec
phi.vector()[:] = phi_vec
Dynamic spherical basis?
e_r_vect = e_r.vector()[:]
e_theta_vect = e_theta.vector()[:]
e_phi_vect = e_phi.vector()[:]
for dof, node in enumerate(vertices(fenics_mesh)):
e_r_vect[dof] = as_vector( (sin(theta_vec[dof])*cos(phi_vec[dof]), sin(theta_vec[dof])*sin(phi_vec[dof]), cos(theta_vec[dof])) ) → setting an array element with a sequence’ ValueError
e_theta_vect[dof] = as_vector( (cos(theta_vec[dof])*cos(phi_vec[dof]), cos(theta_vec[dof])*sin(phi_vec[dof]), sin(theta_vec[dof])) )
e_phi_vect[dof] = as_vector( (-sin(phi_vec[dof]), cos(phi_vec[dof]), 0.) )
e_r.vector()[:] = e_r_vect
e_theta.vector()[:] = e_theta_vect
e_phi.vector()[:] = e_phi_vect
Many thanks!