# Gradient of spatial coordinate wrt element basis?

**URL:** <https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531>\
**Category:** Uncategorized\
**Created:** [February 26, 2020, 4:19pm UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531 "2020-02-26T16:19:26Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![niewiarowski](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/niewiarowski/32/620_2.png) [@niewiarowski](https://fenicsproject.discourse.group/u/niewiarowski)\
**Post date:** [February 26, 2020, 4:19pm UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531/1 "2020-02-26T16:19:26Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [February 26, 2020, 4:21pm UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531/2 "2020-02-26T16:21:18Z")

</div>

Hello,

```auto
from ufl import Jacobian

```

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

---

<div class="post-metadata">

**Author:** ![niewiarowski](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/niewiarowski/32/620_2.png) [@niewiarowski](https://fenicsproject.discourse.group/u/niewiarowski)\
**Post date:** [February 26, 2020, 4:42pm UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531/3 "2020-02-26T16:42:33Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![niewiarowski](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/niewiarowski/32/620_2.png) [@niewiarowski](https://fenicsproject.discourse.group/u/niewiarowski)\
**Post date:** [February 26, 2020, 5:11pm UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531/4 "2020-02-26T17:11:14Z")

</div>

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
```

---

<div class="post-metadata">

**Author:** ![niewiarowski](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/niewiarowski/32/620_2.png) [@niewiarowski](https://fenicsproject.discourse.group/u/niewiarowski)\
**Post date:** [February 27, 2020, 12:07am UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531/5 "2020-02-27T00:07:06Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Thomas2](https://avatars.discourse-cdn.com/v4/letter/t/aeb1de/32.png) [@Thomas2](https://fenicsproject.discourse.group/u/Thomas2)\
**Post date:** [April 1, 2020, 5:20pm UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531/6 "2020-04-01T17:20:29Z")

</div>

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

```auto
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

---

<div class="post-metadata">

**Author:** ![niewiarowski](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/niewiarowski/32/620_2.png) [@niewiarowski](https://fenicsproject.discourse.group/u/niewiarowski)\
**Post date:** [April 14, 2020, 12:27am UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531/7 "2020-04-14T00:27:59Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Thomas2](https://avatars.discourse-cdn.com/v4/letter/t/aeb1de/32.png) [@Thomas2](https://fenicsproject.discourse.group/u/Thomas2)\
**Post date:** [April 15, 2020, 10:05am UTC](https://fenicsproject.discourse.group/t/gradient-of-spatial-coordinate-wrt-element-basis/2531/8 "2020-04-15T10:05:35Z")

</div>

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

```auto
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?