Normal in dolfinx

Karlosjusus · April 21, 2023, 3:36pm

Hi everyone!

I’m trying to project the normal vector

n = FacetNormal(domain)

into a Vector space

v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)

I tried to do somethig like How to plot normal unit vector of faces in a 2D mesh? , but in FEniCSX I dont get it…

Anyone have a idea?

Thank you!

dokken · April 21, 2023, 4:37pm

See: dolfinx_mpc/mpc_utils.py at main · jorgensd/dolfinx_mpc · GitHub

Karlosjusus · April 24, 2023, 5:28pm

I tried to do, but I get the next error.

ERROR:UFL:Invalid subdomain_id [4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 1 4 3 3 3
 1 3 3 3 4 3 4 4 3 4 4 3 4 3 4 4 4 3 4 4 3 4 4 3 4 3 4 3 4 4 3 4 4 4 4 3 4
 4 4 3 4 3 4 4 3 4 4 3 4 4 3 4 4 4 4 3 4 4 3 4 4 3 4 4 3 4 3 4 4 4 4 3 4 3
 3 4 3 4 5 5 5 5 3 5 5 4 3 5 5 5 5 3 4 4 4 5 5 3 5 5 4 4 5 5 4 3 4 5 5 4 3
 5 5 5 5 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 4 4 5 5 5 5
 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 5 5 5 5 5 5
 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
 4 4 5 3 4 5 5 5 5 5 5 5 5 5 5 3 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 5 5 4
 5 5 5 5 5 5 4 5 5 5 5 5 5 5 5 5 5 4 5 5 5 5 5 5 3 4 4 5 5 5 5 5 5 5 5 5 5
 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 5 5 5 5 5 5 5 5 5 5 4 5 5 5 5 5 3 4 4
 5 5 5 5 5 5 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 4 4 5 5 5 5 5 5 5 5 5 5 5
 5 4 4 4 5 5 5 5 5 5 5 5 4 5 5 5 5 5 5 5 5 3 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5
 5 5 5 4 5 5 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 4 4 5
 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 5 5 5 5 5 5 5
 5 5 5 5 3 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 5 5 3 4 5 5 5 5 5 5 5 5 5 5
 4 5 5 5 5 5 5 5 5 5 3 5 5 5 4 5 5 5 3 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 4
 5 5 5 5 5 5 4 5 5 5 5 3 5 5 5 5 5 5 4 5 5 5 5 5 5 5 5 5 3 4 5 5 5 5 5 5 5
 5 5 5 5 5 5 5 5 5 5 3 4 5 5 4 4 4 3 4 3 4 4 3 4 4 3 4 3 4 4 3 4 4 4 4 3 4
 4 4 3 4 4 4 3 4 2 2 2 2 2 2 2 2 2 3 2 4 2 2 2 4 2 4 2 4 2 4 4 4 4 2 4 2].

My code. I tried facet_tag.inidices() too:

import dolfinx
from mpi4py import MPI
import numpy as np
import math
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedElement)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, 
                               create_vector, create_matrix, set_bc)
import matplotlib.pyplot as plt

from dolfinx.fem import (Constant, Function, FunctionSpace, 
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)


##########################################################################################################3

    # Read the mesh

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh_cinta.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh()

    # Definition of Spaces. (dim and polinomial type)

v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)



########################################################################################################

    # Definition of subdomains
    
x_up  = -40 
x_down = 60
y_top = 20
y_floor = 0 
x_left = 0
x_right = 10
y_cavity = -5

n = FacetNormal(domain)

def upstream(x):
    return np.isclose(x[0], x_up)

def downstream(x):
    return np.isclose(x[0], x_down)

def top(x):
    return np.isclose(x[1], y_top)

def cavity_y(x):
    return np.logical_or(np.logical_and(np.logical_or(np.isclose(x[1], y_floor), np.isclose(x[1], y_cavity)), np.logical_or(np.less(x[0], x_left) , np.greater(x[0], x_right))),np.isclose(x[1], y_cavity))

def cavity_x(x):
    return np.logical_and(np.logical_or(np.isclose(x[0], x_left), np.isclose(x[0], x_right)),  np.less(x[1], 0.1))

    # Whole cavity
    
def cavity(x):
    return np.logical_or(cavity_x(x) , cavity_y(x))

def cinta(x):
    return np.logical_and( np.logical_and(np.greater(x[0], x_left) , np.less(x[0], x_right)) , np.logical_and(np.greater(x[1],y_cavity) , np.less(x[1],5)) )



#############################################################################################################

    #Define normal and tangent like functions

def facet_normal_approximation(V, mt: dolfinx.mesh.meshtags, mt_id: int, tangent=False, jit_options: dict = {},
                               form_compiler_options: dict = {}):

    timer = dolfinx.common.Timer("~MPC: Facet normal projection")
    comm = V.mesh.comm
    n = ufl.FacetNormal(V.mesh)
    nh = Function(V)
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    ds = ufl.ds(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id)
    if tangent:
        if V.mesh.geometry.dim == 1:
            raise ValueError("Tangent not defined for 1D problem")
        elif V.mesh.geometry.dim == 2:
            a = ufl.inner(u, v) * ds
            L = ufl.inner(ufl.as_vector([-n[1], n[0]]), v) * ds
        else:
            def tangential_proj(u, n):
                """
                See for instance:
                https://link.springer.com/content/pdf/10.1023/A:1022235512626.pdf
                """
                return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u
            c = dolfinx.fem.Constant(V.mesh, [1, 1, 1])
            a = ufl.inner(u, v) * ds
            L = ufl.inner(tangential_proj(c, n), v) * ds
    else:
        a = (ufl.inner(u, v) * ds)
        L = ufl.inner(n, v) * ds

    # Find all dofs that are not boundary dofs
    imap = V.dofmap.index_map
    all_blocks = np.arange(imap.size_local, dtype=np.int32)
    top_blocks = dolfinx.fem.locate_dofs_topological(V, V.mesh.topology.dim - 1, mt.find(mt_id))
    deac_blocks = all_blocks[np.isin(all_blocks, top_blocks, invert=True)]

    # Note there should be a better way to do this
    # Create sparsity pattern only for constraint + bc
    bilinear_form = _fem.form(a, jit_options=jit_options,
                              form_compiler_options=form_compiler_options)
    pattern = dolfinx.fem.create_sparsity_pattern(bilinear_form)
    pattern.insert_diagonal(deac_blocks)
    pattern.assemble()
    u_0 = _fem.Function(V)
    u_0.vector.set(0)

    bc_deac = dolfinx.fem.dirichletbc(u_0, deac_blocks)
    A = dolfinx.cpp.la.petsc.create_matrix(comm, pattern)
    A.zeroEntries()

    # Assemble the matrix with all entries
    form_coeffs = dolfinx.cpp.fem.pack_coefficients(bilinear_form)
    form_consts = dolfinx.cpp.fem.pack_constants(bilinear_form)
    _cpp.fem.petsc.assemble_matrix(A, bilinear_form, form_consts, form_coeffs, [bc_deac])
    if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
        A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
        A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)
        dolfinx.cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [bc_deac], 1.0)
    A.assemble()
    linear_form = dolfinx.fem.form(L, jit_options=jit_options,
                            form_compiler_options=form_compiler_options)
    b = dolfinx.fem.petsc.assemble_vector(linear_form)

    dolfinx.fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    dolfinx.fem.petsc.set_bc(b, [bc_deac])

    # Solve Linear problem
    solver = PETSc.KSP().create(MPI.COMM_WORLD)
    solver.setType("cg")
    solver.rtol = 1e-8
    solver.setOperators(A)
    solver.solve(b, nh.vector)
    nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    timer.stop()
    return nh

up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
cavity_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, cavity)
top_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, top)
cinta_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, cinta)  

boundaries = [(1,up_facets),(2,down_facets),(3,top_facets),(4,cavity_facets),(5,cinta_facets)]

facet_indices, facet_markers = [], []
for (marker, face ) in boundaries:
    facet_indices.append(face)
    facet_markers.append(np.full_like(face, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

nh = facet_normal_approximation(V, facet_tag, facet_tag.values, tangent=False)

The mesh

dokken · April 24, 2023, 6:09pm

The third argument here should not be a list, But either:

A single integer
A tuple of unique integers. i.e. tuple(np.unique(mt.values))

Karlosjusus · April 24, 2023, 6:28pm

Thank you so much again, @dokken

Topic		Replies	Views
How to plot normal unit vector of faces in a 2D mesh?	7	3932	August 10, 2020
How to compute the normal fector field of a surface embedded in 3D, in FEniCSx? dolfinx	4	303	January 20, 2023
Manually compute normal vector on inner boundary/subdomain	15	1747	November 30, 2019
Normal vector in manifold General	4	147	April 10, 2024
Extract normal vector data at boundary in dolfinx dolfinx	6	330	March 24, 2024

Normal in dolfinx

Related topics