How to compute tangent vectors on facets (2D surface parametrization)

[deshik]

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

[hherlyng]

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

[deshik]

The facetjacobian function in UFL helped me!
from ufl.classes import FacetJacobian, FacetJacobianInverse

Then i use them in residual of my variational problem.

[hussleJ6]

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;

    # 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;

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

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

[hussleJ6]

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;

    # 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;

    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;

Traceback (most recent call last):
  File "/home/username/normals_calc.py", line 240, in <module>
    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.

[hussleJ6]

Ok, I looked somewhat further and bumped into a post with the exact same issue;

Cannot write function to XDMF with DG elements

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

Hi all, This is a simpler version of the problem stated in Imposing a discontinuity in the solution, on an internal interface 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;

from dolfinx import default_scalar_type

but it is rather imported with;

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
For what reason am I encountering this error in that case;

Traceback (most recent call last):
  File "/home/username/normals.py", line 223, in <module>
    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

    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?

[dokken]

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 with XDMFFile there was an implicit interpolation step into continuous Lagrange prior to storing it.

[hherlyng]

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 you should be free of NaN values.

Best,
Halvor