Hello everyone.
I am trying to project several functions, obtained from external data. When I project some functions, I get that they are 0, but the original functions are not. In fact there is one function, which does not give zero, but I don’t understand why?
Chi_expr is equal to zero but S_expr no.
Any idea?
In this folder is the mesh and the necessary data for the code to compile.
my version is 0.7.3
import dolfinx
import gmsh
from mpi4py import MPI
from petsc4py import PETSc
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)
from dolfinx.mesh import create_unit_square
from petsc4py.PETSc import ScalarType
from dolfinx import fem
from dolfinx import geometry
from dolfinx.io import gmshio
#import meshio
import dolfinx.cpp as _cpp
import dolfinx.fem as _fem
from dolfinx.nls.petsc import NewtonSolver
print(f"DOLFINx version: {dolfinx.__version__}")
domain, cell_tags, ft = gmshio.read_from_msh("Mesh/cylinder.msh", MPI.COMM_WORLD, 0, gdim=2)
ft.name = "Facet markers"
#Encontrar los id segun las coordenadas
tdim = domain.topology.dim
tree = dolfinx.geometry.bb_tree(domain, tdim)
num_entities_local = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts
entities = np.arange(num_entities_local, dtype = np.int32)
mid_tree = dolfinx.geometry.create_midpoint_tree(domain, tdim,entities)
mesh_geom = domain.geometry.x
tol = 1e-8
def closest_point(mesh, point,tol):
points = point
while True:
try:
entity = dolfinx.geometry.compute_closest_entity(tree, mid_tree, mesh, points)
break
except RuntimeError:
print(points)
for j in range(len(mesh.geometry.x)):
if (np.abs(mesh_geom[j][0]-points[0])< tol and np.abs(mesh_geom[j][1]-points[1])< tol):
return points, j
geom_dof = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object,tdim,entity, False)
#geom_dof = dolfinx.cpp.mesh.compute_incident_entities(mesh,[entity], 2,1)
#print('geom',geom_dof)
mesh_nodes = mesh_geom[geom_dof][0]
#print('mesh_nodes', mesh_nodes)
dis = dolfinx.geometry.compute_distance_gjk(points, mesh_nodes)
index = -100
for i in range(len(mesh_nodes)):
#print(mesh_nodes[i])
if (np.abs(mesh_nodes[i][0]-points[0]+ dis[0]) < tol and np.abs(mesh_nodes[i][1]-points[1]+ dis[1]) < tol):
index = i
if (index == -100):
return point, index
else:
return points[0] - dis[0], points[1] - dis[1] , geom_dof[0][index]
########################################################################################################
# Definition of Spaces. (dim and polinomial type)
v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
v_cg2 = VectorElement("CG", domain.ufl_cell(), 1)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
s_cg2 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg2)
Q = dolfinx.fem.FunctionSpace(domain, s_cg2)
Z = MixedElement([v_cg1, s_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)
V1, _ = W.sub(0).collapse()
Q1, _ = W.sub(1).collapse()
dataname = 'data.txt'
vx_np = np.loadtxt(dataname)
vx_aux = vx_np[:,0:2]
b = []
indices = []
ind = int(0)
for i in range(len(vx_aux)):
a,y,j = closest_point(domain,vx_aux[i],0.0001)
indices.append(int(j))
print(len(vx_aux))
data_np = np.loadtxt(dataname)
u = Function(V1)
nu_t1 = Function(Q)
d = Function(Q)
u_aux = Function(V)
# AnsysFluent to FEniCSX
###################################################################################################################
# 0 = x_coordinate 1 = y_coordinate 2 = nu turbulenta 3 = nu effective 4 = x_velocity 5 = y_velocity 6 = density 7 = pressure #
###################################################################################################################
for i in range(len(indices)):
nu_t1.x.array[indices[i]] = data_np[i,2]
u_aux.x.array[2*indices[i]+1]=data_np[i,5]
d.x.array[indices[i]] = data_np[i,8]
u_aux.x.array[2*indices[i]]=data_np[i,4]
u.interpolate(u_aux)
nu_t = Function(Q1)
nu_t.interpolate(nu_t1)
# Constantes del modelo
Cv1 = Constant(domain, PETSc.ScalarType(7.1))
nu_p = Constant(domain, PETSc.ScalarType(9.4e-5))
# Lista de valores para nu_t
nu_t_list = nu_t.x.array[:] # Aquí colocas tus valores para nu_t
nu = 9.4e-5 # Valor de la viscosidad cinemática molecular
C_v1 = 7.1
# Función para resolver el polinomio para cada nu_t
def calcular_tilde_nu(nu_t, nu, C_v1):
# Coeficientes del polinomio de cuarto grado
coeffs = [1, -nu_t, 0, 0, -nu_t * (C_v1**3)*nu**3]
# Encontrar las raíces del polinomio
roots = np.roots(coeffs)
# Seleccionar la raíz real positiva
tilde_nu = np.real(roots[roots > 0])
# Selecciona la raíz más pequeña real positiva
return tilde_nu[0] if len(tilde_nu) > 0 else None
# Aplicar la función a cada valor de la lista nu_t
tilde_nu_list = [calcular_tilde_nu(nu_t, nu, C_v1) for nu_t in nu_t_list]
nu_tilde = Function(Q)
nu_tilde.x.array[:] = tilde_nu_list
grad_u = ufl.grad(u)
Omega = 0.5 * (grad_u - ufl.transpose(grad_u))
S_expr = ufl.sqrt(2*ufl.inner(Omega, Omega))
chi_expr = nu_tilde / nu
def project_expression(expr, space):
v = ufl.TestFunction(space)
u_trial = ufl.TrialFunction(space)
a = ufl.inner(u_trial, v) * ufl.dx
L = ufl.inner(expr, v) * ufl.dx
problem = fem.petsc.LinearProblem(a, L)
return problem.solve()
# Proyectamos cada una de las funciones auxiliares en el espacio S_space:
chi_proj = project_expression(chi_expr, Q)
S_proj = project_expression(S_expr, Q)
def write_function_to_file(func, filename,):
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(func)
write_function_to_file(chi_proj, "Plot/chi.xdmf")
write_function_to_file(S_proj, "Plot/S.xdmf")
