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!
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
The third argument here should not be a list, But either:
tuple(np.unique(mt.values))Thank you so much again, @dokken