Gradient of spatial coordinate wrt element basis?

1

Hi all,
Given a manifold mesh (topological dimension=2, geometric dimension=3), for example the surface mesh of the UnitCubeMesh,

from dolfin import *
mesh = UnitCubeMesh(2,2,2)
srf = BoundaryMesh(mesh, "exterior")
srf.topology().dim()  # is 2
srf.geometric_dimension()  # is 3
x = SpatialCoordinate(srf)

how can I obtain the UFL representation of the element tangent basis (the jacobian components of the mapping from the reference domain)? Explicitly, I would like: \partial {x} / \partial \xi^1 and \partial x / \partial \xi^2 , where \xi is the parent domain coordinate (or alternatively, the shape function derivatives wrt to the reference domain coordinates \partial N_I /\partial \xi^I. And would this be possible for triangular elements in addition to quads? Thanks in advance!

2

Hello,

from ufl import Jacobian

see also https://fenics.readthedocs.io/projects/ufl/en/latest/api-doc/ufl.html#ufl.Jacobian

3

Hi @bleyerj, thanks for the quick response - is it possible to create \partial {x} / \partial \xi^1 and \partial x / \partial \xi^2 separately, or extract them using Index? I tried below to plot the \xi^1 tangent basis but it doesn’t seem to work as expected.

from dolfin import *
mesh =UnitCubeMesh(2,2,2)
srf = BoundaryMesh(mesh, "exterior")
x = SpatialCoordinate(srf)
from ufl import Jacobian, Index
J = Jacobian(srf)
i = Index()
G1 = project(J[i,0], VectorFunctionSpace(srf, 'DG', 1, dim=3))

from vtkplotter.dolfin import plot
plot(G1)
4

The MWE below “works”, but doesn’t give me the result I expect (G1xG2 = non normalized outward facing normals). How can I orient the cells such that the calculated normal is outward facing? Thanks for your help

from dolfin import *
mesh =UnitCubeMesh(3,3,3)
srf = BoundaryMesh(mesh, "exterior")
global_normal = Expression(("x[0]", "x[1]", "x[2]"), degree = 1)
srf.init_cell_orientations(global_normal)

x = SpatialCoordinate(srf)
from ufl import Jacobian
J = Jacobian(srf)
G1 = J[:,0]
G2 = J[:,1]
N = cross(G1, G2)
Nh = project(N, VectorFunctionSpace(srf, 'DG', degree=2, dim=3))
File('Nh.pvd') << Nh
5

Marking this as solved, thanks. Will post my followup separately

6

Hello @niewiarowski
I am also interested in creating the Jacobian directly. In your last post mentioned a “followup” , I found
"Consistent non-normalized normal field for immersed ", which addresses the direction of the normal . Is what you meant by “followup” ?

But perhaps you have found a way to create the Cell Jacobian directly? Please let me know if you have.

Some background about why this would be useful:
For our application we have an \mathbb R^n manifold immersed in \mathbb R^m (in this case n=2 and m=3) , the standard Jacobian function in UFL returns
J_{ij}=\frac{\partial x_i}{\partial X_j} for i=1....m and j=1...n
but for our application we need (or at least we think we need…)
J_{ij}=\frac{\partial x_i}{\partial X_j} for i=1....m and j=1...m
where we manually create X_m from the normal of the manifold

so far we have added an extra function to ffc/backends/ufc/ufc_geometry.h in order to calculate this Jacobian, our C++ code looks like this

inline void crossv1v2(double ret[UFC_GDIM_3],const double v1[3],const double v2[3]){
    for (int i=0;i<3;i++)ret[i]=v1[(i+1)%3]*v2[(i+2)%3]-v1[(i+2)%3]*v2[(i+1)%3];
};

inline void compute_jacobian_triangle_3d_Cosserat(double J[UFC_GDIM_3*UFC_GDIM_3], const double coordinate_dofs[9]){
  double cellNormal[UFC_GDIM_3],v1[UFC_GDIM_3],v2[UFC_GDIM_3];
  for (int i=0;i<3;i++){
      v1[i]=coordinate_dofs[i+3]-coordinate_dofs[i];
      v2[i]=coordinate_dofs[i+3+3]-coordinate_dofs[i+3];
  };
  crossv1v2(cellNormal,v1,v2);
  J[0] = coordinate_dofs[3]  - coordinate_dofs[0];
.....similar lines removed for brevity
  J[8] = cellNormal[2] - coordinate_dofs[2];
}

but we aren’t sure how to link this C++ function to make it available as ufl a python function

(Actually our “manifolds” aren’t really manifolds because they are only “piece-wise differentiable” rather than differentiable everywhere but I don’t think thats relevant here)

Regards
Thomas

7

Hi @Thomas2

I just saw your response to my other question - all this is for a side project I was working on and I haven’t looked at it in a few weeks, but I’m taking a second swing at it. Have you had any luck? By the way, it should be possible to create python bindings for your custom function

8

Hi Alexander

we gave up on the C++ route , but we have a “sort of working” UFL version - we are attempting to transform the global unknown in \mathbb R^3 into the tangent space + normal -( so just a rotated version of \mathbb R^3 )
We create the basis vectors in ufl, then iterate through each dimension and use the ufl.dot product with the global unknown-

The code is a bit like this

e1=ufl.geometry.CellEdgeVectors(domain)[0,...]
e3=ufl.CellNormal(domain)
e2=ufl.cross(e1,e3); 
basis=[-e1,-e2,-e3] 

for iv ,v in enumerate(basis):
    basis[iv]/=ufl.sqrt(ufl.dot(v,v))

eltDeformationRotation=[ufl.dot(globalDeformation,basis[ii])  for ii in range(3)]+[ufl.dot(globalRotation,basis[ii]) for ii in range(3)]
eltDeformation=ufl.as_tensor(eltDeformationRotation[:3])
global_grad_eltDeformation=ufl.grad(eltDeformation)
grad_deformation_in_elt_manifold=ufl.as_tensor([[ufl.dot(global_grad_eltDeformation[ii,...],basis[jj]) for jj in range(3)] for ii in range(3)]  )

We also tried creating a metric with ufl.indices , it did work but for some reason the FFC/Dijitso part took ages when we used that approach - probably some error in our code.

Hope this is helpful , please don’t place too much weight on this guidance, I am not very familiar with fenics or the maths…

One thing I have found useful is David Kamensky’s tIGAr - it helped me to work out how to get UFL to do these kind of operations - probably that would be useful to your work as well?

1 Like