# How to compute tangent vectors on facets (2D surface parametrization)
**URL:** https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761
**Category:** General
**Created:** [November 15, 2023, 7:59pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761 "2023-11-15T19:59:27Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@deshik](https://fenicsproject.discourse.group/u/deshik)
#### Post date: [November 15, 2023, 7:59pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761/1 "2023-11-15T19:59:27Z")
Hi,
I have a 3d mesh with two domains and an interface. I would like to do surface parametrization on the interface. For this I need to compute tangent vectors on the interface facets.
That is,$$ \partial x / ‘\partial X’ $$
x = physical coordinates,
X = reference coordinates.
I think I can compute this just using the [Jacobian](https://fenics.readthedocs.io/projects/ufl/en/2019.1.0/api-doc/ufl.html#ufl.classes.Jacobian) function for 3D tetrahedra, but i couldn’t calculate it for facet.
Can someone help me how i can do this?
Just FYI, I have tagged the domains, boundaries and interface using GMSH, and I can successfully read them in fenics.
Thanks!
---
### Author:  [@hherlyng](https://fenicsproject.discourse.group/u/hherlyng)
#### Post date: [November 17, 2023, 1:15pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761/3 "2023-11-17T13:15:07Z")
Hi deshik,
Since there is no detailed specification of what programming language/package you are using to attempt this I went ahead and assumed that an implementation based on dolfinx is fine for you.
In [this gist](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f) you can find a code for approximating facet vectors on all facets or a subset of facets on a dolfinx mesh. The function has a flag `tangent`. When this flag is set to `False` the facet normal vectors will be approximated, whereas the facet tangent vectors are approximated when the flag is set to `True`.
The code also has a flag `interior`. Set this to `True` if the facet vectors you are approximating are on interior facets, such as in your case for the interface between two domains. Set the flag to `False` if you want to approximate facet vectors on the boundary of the mesh.
The facet tags of the facets you want to approximate the facet vectors on are passed in as the argument `mt`, with the corresponding tag value `mt_id` of those meshtags.
Note that the approximated vectors are defined in a DG1 finite element space. Resultingly, visualizing output files written in the script will display several vectors in each vertex since there are facet vectors defined on facets of all neighboring cells of the vertex. In practice however when e.g. integrating over a facet that facet’s normal/tangent vector is uniquely determined.
Hope this is of help!
Cheers,
Halvor
---
### Author:  [@deshik](https://fenicsproject.discourse.group/u/deshik)
#### Post date: [January 17, 2024, 7:42pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761/4 "2024-01-17T19:42:21Z")
The facetjacobian function in UFL helped me!
from ufl.classes import FacetJacobian, FacetJacobianInverse
Then i use them in residual of my variational problem.
---
### Author:  [@hussleJ6](https://fenicsproject.discourse.group/u/hussleJ6)
#### Post date: [July 13, 2024, 11:42am UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761/5 "2024-07-13T11:42:06Z")
Dear @hherlyng, I had a question with regards to this code which calculates the normals to facets. As soon as the normals.bp file is written, I try to visualize and extract the normals from this file. When the normals.bp file is opened in ParaView, 3D glyphs only show vectors in the x-direction. For my mesh, I obtain the exact same thing. So to visualize the normal data, I have tried to write the normals into a .txt file. Furthermore I want to obtain the normalized normal vectors on every single node, to multiply those with a scalar value. Apparently it obtains a lot of nan values. In what way could the normals be extracted at every single node (or at least on every facet)? Why do the nan values occur and should those nan values be set to zero for example? I would expect a format like 1 1 1. Right now it is in such a format, if written wit this code;
```auto
# Assume nh is a dolfinx.fem.Function
nh_values = nh.vector[:]
# Define the filename for the TXT output
txt_filename = "normal_values_only.txt"
# Write data to TXT file
np.savetxt(txt_filename, nh_values, header='NH', comments='', fmt='%f')
```
Format;
```auto
-0.316228
0.000000
-0.948683
0.707107
0.000000
-0.707107
0.707107
0.000000
-0.707107
0.948683
0.000000
0.316228
-0.316228
-0.948683
0.000000
0.707107
-0.707107
0.000000
0.948683
0.316228
0.000000
0.707107
-0.707107
0.000000
0.000000
0.000000
-1.000000
0.000000
0.000000
-1.000000
0.000000
0.000000
-1.000000
nan
nan
nan
0.000000
-1.000000
0.000000
0.000000
-1.000000
0.000000
nan
nan
nan
0.000000
-1.000000
0.000000
nan
nan
nan
1.000000
0.000000
0.000000
1.000000
0.000000
0.000000
1.000000
0.000000
0.000000
nan
nan
nan
nan
nan
nan
```
Those were set;
```auto
all_facets = True
tangent_flag = False
dim = 3
facet_tags = None
ft_id = None
```
Is it possible to obtain the normal vectors at every single node and write it to a .txt file or a .xdmf/.h5 file? I would like to hear from you.
This is how it looks in ParaView if the .bp file is loaded and viewed with 3D\_Glyphs
---
### Author:  [@hussleJ6](https://fenicsproject.discourse.group/u/hussleJ6)
#### Post date: [July 16, 2024, 9:40am UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761/6 "2024-07-16T09:40:13Z")
Hi there, I had a question on saving the normal vector solution nh in another format. In the code which was provided, the solution has been written to a .bp format. I would prefer to write it into the format of a .xdmf/.h5 file, but there is a constant discrepancy between both files. Probably this is due to the fact that a “Continuous Lagrange” family was chosen to write the .xdmf file, in contrast to the .bp file, were a “Discontinuous Lagrange” family was chosen. In the images down below is shown the differences for both formats. This was the original code;
```auto
# Create a DG1 space for the facet vectors to be approximated.
DG1 = ufl.VectorElement(family='Discontinuous Lagrange', cell=mesh.ufl_cell(), degree=1)
space = dfx.fem.FunctionSpace(mesh=mesh, element=DG1)
# Compute the facet vector approximation.
nh = facet_vector_approximation(V=space, mt=facet_tags, mt_id=ft_id, interior=interior, tangent=tangent_flag)
# Write the results to file
filename = 'tangents.bp' if tangent_flag else 'normals.bp'
with dfx.io.VTXWriter(mesh.comm, filename, [nh], engine='BP4') as vtx:
vtx.write(0)
```
To write it to a .xdmf file, the following was added;
```auto
xdmf_filename = "normals.xdmf"
# Create a continuous vector space for visualization
CG1 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), degree=1)
V_vis = dfx.fem.FunctionSpace(mesh, CG1)
# Interpolate discontinuous data to the continuous space
nh_vis = dfx.fem.Function(V_vis)
nh_vis.interpolate(nh)
# Now write to XDMF
with dolfinx.io.XDMFFile(mesh.comm, xdmf_filename, "w") as xdmf_file:
xdmf_file.write_mesh(mesh)
xdmf_file.write_function(nh_vis)
```
So probably the difference in Lagrange family type is causing this to happen, but if “Discontinuous Lagrange” was chosen for the last as well, this error occurs;
```auto
Traceback (most recent call last):
File "/home/username/normals_calc.py", line 240, in
xdmf_file.write_function(nh_vis)
File "/home/username/miniconda3/envs/env_name/lib/python3.10/site-packages/dolfinx/io/utils.py", line 235, in write_function
super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)
RuntimeError: Function and Mesh dof layouts do not match. Maybe the Function needs to be interpolated?
```
Does anyone know how to resolve this issue? Thanks for any help already.
---
### Author:  [@hussleJ6](https://fenicsproject.discourse.group/u/hussleJ6)
#### Post date: [July 17, 2024, 1:10pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761/7 "2024-07-17T13:10:14Z")
Ok, I looked somewhat further and bumped into a post with the exact same issue;
> [@Cannot write function to XDMF with DG elements](https://fenicsproject.discourse.group/t/cannot-write-function-to-xdmf-with-dg-elements/12442):
>
> Hi! I’m trying to export a function from a DG function space to XDMF but obtain the following error: File "/home/remidm/festimx/mwe.py", line 14, in \ xdmf\_file.write\_function(u) File "/home/remidm/miniconda3/envs/dolfinx/lib/python3.11/site-packages/dolfinx/io/utils.py", line 235, in write\_function super().write\_function(getattr(u, "\_cpp\_object", u), t, mesh\_xpath) RuntimeError: Function and Mesh dof layouts do not match. Maybe the Function needs to be interpolated? Obtain…
As it seems, the current dolfinx version does not support to write the solution to a .xdmf file when a ‘Discontinuous Lagrange’ space was chosen. Whereas a ‘Continuous Lagrange’ space does not seem to be sufficient in this particular case, rather a possibility would be to switch to write the solutions with ADIOS & make use of the VTXwriter. Despite I’ ve also seen that it has been able to write the solution to a .xdmf file once was chosen to make use of the older dolfinx-0.7.0 version.
> [@Imposing a discontinuity at interface using DG method](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-at-interface-using-dg-method/11780):
>
> Hi all, This is a simpler version of the problem stated in [Imposing a discontinuity in the solution, on an internal interface](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-in-the-solution-on-an-internal-interface/10112) I’d like (in dolfinx) to impose a discontinuity at the interface between two domains. I understood that MeshView would be the best option since we would have CG elements everywhere except at the interface. However it’s not yet implemented in FEniCSx I’d like to give the DG method a go! This is a simple poisson problem and I’d like to enforce the following equation at…
I would still prefer to write the solution to a .xdmf file and (luckily for me) am still working in this old 0.7.0 version.
Despite, a couple of errors are encountered and I was wondering what might be the cause. For example default\_scalar type is not imported anymore in those legacy versions with;
```auto
from dolfinx import default_scalar_type
```
but it is rather imported with;
```auto
from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
```
Probably there are a couple of those changes. What i don’t understand is that if I look in the documentation of this old dolfinx version, it already seems to support for example finalize;
[https://docs.fenicsproject.org/dolfinx/v0.7.0/python/generated/dolfinx.cpp.la.html?highlight=finalize#dolfinx.cpp.la.SparsityPattern.finalize](https://docs.fenicsproject.org/dolfinx/v0.7.0/python/generated/dolfinx.cpp.la.html?highlight=finalize#dolfinx.cpp.la.SparsityPattern.finalize)
For what reason am I encountering this error in that case;
```auto
Traceback (most recent call last):
File "/home/username/normals.py", line 223, in
nh = facet_vector_approximation(V=space, mt=facet_tags, mt_id=ft_id, interior=interior, tangent=tangent_flag)
File "/home/username/normals.py", line 70, in facet_vector_approximation
pattern.finalize()
AttributeError: 'dolfinx.cpp.la.SparsityPattern' object has no attribute 'finalize'
```
for this piece of code
```auto
bilinear_form = dfx.fem.form(a, jit_options=jit_options,
form_compiler_options=form_compiler_options)
pattern = dfx.fem.create_sparsity_pattern(bilinear_form)
pattern.insert_diagonal(deac_blocks)
pattern.finalize()
```
To me this does not seem to make sense. Did anyone experience something similar?
Edit: Changing this to pattern.assemble() & removing every \_cpp\_object made it possible to write the solution to a .xdmf format with a Discontinuous Lagrange space. Still the pattern with the spots remains, like the case when the continuous Lagrange was chosen as before. Does anyone know why this is exactly the case. Why does it leave some spots empty and does not allocate a value to some nodes, although it does when using VTXwriter?
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [August 25, 2024, 6:58pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761/8 "2024-08-25T18:58:56Z")
The XDMFFile format has not written anything other than first order continuous interpolations for a very long time. As stated in: [Cannot write function to XDMF with DG elements - #2 by dokken](https://fenicsproject.discourse.group/t/cannot-write-function-to-xdmf-with-dg-elements/12442/2) with XDMFFile there was an implicit interpolation step into continuous Lagrange prior to storing it.
---
### Author:  [@hherlyng](https://fenicsproject.discourse.group/u/hherlyng)
#### Post date: [August 29, 2024, 3:09pm UTC](https://fenicsproject.discourse.group/t/how-to-compute-tangent-vectors-on-facets-2d-surface-parametrization/12761/9 "2024-08-29T15:09:03Z")
Hi hussleJ6,
sorry for the belated reply over the summer. I suspect the NaN values were due to a missing check as to whether the length of a normal vector in a point was zero, which resulted in division by zero during normalization of the vectors. A check has been added to the code now, so using [this updated version of the code](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f)you should be free of NaN values.
Best,
Halvor