How to increase number of boundary points?

General
1

Hello
I am using to following code for creating the mesh model based on the the code in the Dokken’s FEniCSx tutorial for the flow past a cylinder (Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial ).

############################################
#--------------- Imports ------------------#
############################################

import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import tqdm.autonotebook
import matplotlib as mpl

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_vector, create_matrix, set_bc)
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities

from dolfinx import plot

from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)

import dolfinx
import viskex
import time
import pylab as pl
from IPython import display
import matplotlib.tri as tri

mesh_comm = MPI.COMM_WORLD

############################################
#--------------- Mesh creation ------------#
############################################
model_rank = 0
gmsh.initialize()

# L = 2.2
# H = 0.41
# c_x = c_y = 0.2
# r = 0.05
L = 3e-4
H = 5e-5
c_x , c_y = (5e-5, 25e-6)
r = 25e-7
gdim = 2

## Creation of the geometry entities (a cylinder within a rectangle)
if mesh_comm.rank == model_rank:
    rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

## The next step is to substract the obstacle from the channel, such that we do not mesh the interior of the circle
if mesh_comm.rank == model_rank:
    fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
    gmsh.model.occ.synchronize()

## To get GMSH to mesh the fluid, we add a physical volume marker
fluid_marker = 1
if mesh_comm.rank == model_rank:
    volumes = gmsh.model.getEntities(dim=gdim)
    assert (len(volumes) == 1)
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
    gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

## To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side)
## with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by computing the center of mass for each geometrical
## entity.

inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, H / 2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H / 2, 0]):
            outflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
            walls.append(boundary[1])
        else:
            obstacle.append(boundary[1])
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

res_min = 5e-6 #for micro
# res_min =r/1.5
distance_field = gmsh.model.mesh.field.add("Distance")
edges = gmsh.model.getBoundary(volumes, oriented=False)
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", [e[1] for e in edges])
threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax",res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", res_min)
min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)

### GENERATING THE MESH
## We are now ready to generate the mesh. However, we have to decide if our mesh should
## consist of triangles or quadrilaterals. Here we use second order quadrilateral elements.

if mesh_comm.rank == model_rank:
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 1)
    gmsh.option.setNumber("Mesh.RecombineAll", 1)
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")

#=============== Loading mesh and boundary markers============================
## As we have generated the mesh, we now need to load the mesh and corresponding
## facet markers into DOLFINx. See https://jsdokken.com/src/tutorial_gmsh.html for explanations
mesh, c, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
mesh.name = "Cylinder"
c.name = "Cell markers"
ft.name = "Facet markers"

how I could increase the number of boundary points i.e inlet_marker, outlet_marker, wall_marker, and obstacle_marker without decreasing mesh size?

2

How would you be able to not change the mesh size, but introduce more points at the boundary? With more points you have to introduce additional cells to the mesh.

You can of course modify the distance fields specified in:

to change the mesh size distance from edges with given resolution etc.

3

I do not know how I should change the code to get more points currently for inlet_marker, outlet_marker, wall_marker, and obstacle_marker are 21, 21, 242, and 16 points. I want to increase number of points to 128, 128, 1024 and 256.
Thanks

4

This is a GMSH question if you want to do it at the time of meshing.

Then you can adjust the gmsh fields (as shown in the example above and in Gmsh 4.12.2

Alternatively, you could do the mesh refinement in dolfinx, using refine_plaza (which can transfer mesh tags to the refined mesh). For that refinement, you use dolfinx.mesh.locate_entities_boundary to mark the edges of the mesh you want to refine.
An example of such refinement can be found at: dolfinx/python/test/unit/mesh/test_refinement.py at 52348199a9b95c0e38d3e4604505bc68bcc275eb · FEniCS/dolfinx · GitHub

5

I got the following error by using the code:

mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
boundary_facets = ft.find(obstacle_marker)
boundary_vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, boundary_facets, 1, 0)
fine_mesh, parent_cell, parent_facet = dolfinx.mesh.refine_plaza(mesh, boundary_vertices, False, dolfinx.mesh.RefinementOption.parent_cell_and_facet)

RuntimeError: Cell type not supported

6

Right, dolfinx does not support refinement pg quadrilaterals, as you end up with uniform strips of refinement. I would then use Gmsh to tweak the resolution at a bondary

7

Yeah but using Gmesh would increase number of points in all areas. for instance in the code below I just want to increase the points in the obstacle region so the LcMin is decreased. By lowering the mesh size the points in all areas increased but I want only increasing in the obstcle area.

res_min = 5e-6

gmsh.model.mesh.field.add("Distance",1)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", obstacle)
gmsh.model.mesh.field.add("Threshold",2)
gmsh.model.mesh.field.setNumber(2, "InField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", r/3)
gmsh.model.mesh.field.setNumber(2, "LcMax", 0.25 * H)
gmsh.model.mesh.field.setNumber(2, "DistMin", r)
gmsh.model.mesh.field.setNumber(2, "DistMax", 4 * H)

gmsh.model.mesh.field.add("Distance",3)
gmsh.model.mesh.field.setNumbers(3, "EdgesList", outflow)
gmsh.model.mesh.field.add("Threshold",4)
gmsh.model.mesh.field.setNumber(4, "InField", 3)
gmsh.model.mesh.field.setNumber(4, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(4, "LcMax", 0.25 * H)
gmsh.model.mesh.field.setNumber(4, "DistMin", r)
gmsh.model.mesh.field.setNumber(4, "DistMax", 4 * H)

gmsh.model.mesh.field.add("Distance",5)
gmsh.model.mesh.field.setNumbers(5, "EdgesList", inflow)
gmsh.model.mesh.field.add("Threshold",6)
gmsh.model.mesh.field.setNumber(6, "InField", 5)
gmsh.model.mesh.field.setNumber(6, "LcMin", r)
gmsh.model.mesh.field.setNumber(6, "LcMax", 0.25 * H)
gmsh.model.mesh.field.setNumber(6, "DistMin", r)
gmsh.model.mesh.field.setNumber(6, "DistMax", 4 * H)

gmsh.model.mesh.field.add("Distance",7)
gmsh.model.mesh.field.setNumbers(7, "EdgesList", walls)
gmsh.model.mesh.field.add("Threshold",8)
gmsh.model.mesh.field.setNumber(8, "InField", 7)
gmsh.model.mesh.field.setNumber(8, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(8, "LcMax", 0.25 * H)
gmsh.model.mesh.field.setNumber(8, "DistMin", r)
gmsh.model.mesh.field.setNumber(8, "DistMax", 4 * H)

min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [2, 4, 6, 8])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)
8

This just means that you should change LC max to a bigger value (and probably distmax to a smaller value).

9

I changed the code to:

res_min = 5e-6 #for micro
# res_min =r/1.5
gmsh.model.mesh.field.add("Distance",1)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", obstacle)
gmsh.model.mesh.field.add("Threshold",2)
gmsh.model.mesh.field.setNumber(2, "InField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", 5e-7)
gmsh.model.mesh.field.setNumber(2, "LcMax", 10 * H)
gmsh.model.mesh.field.setNumber(2, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(2, "DistMax",  res_min)

gmsh.model.mesh.field.add("Distance",3)
gmsh.model.mesh.field.setNumbers(3, "EdgesList", outflow)
gmsh.model.mesh.field.add("Threshold",4)
gmsh.model.mesh.field.setNumber(4, "InField", 3)
gmsh.model.mesh.field.setNumber(4, "LcMin", r/3)
gmsh.model.mesh.field.setNumber(4, "LcMax", 7e-6)
gmsh.model.mesh.field.setNumber(4, "DistMin", r)
gmsh.model.mesh.field.setNumber(4, "DistMax", r)

gmsh.model.mesh.field.add("Distance",5)
gmsh.model.mesh.field.setNumbers(5, "EdgesList", inflow)
gmsh.model.mesh.field.add("Threshold",6)
gmsh.model.mesh.field.setNumber(6, "InField", 5)
gmsh.model.mesh.field.setNumber(6, "LcMin", r/3)
gmsh.model.mesh.field.setNumber(6, "LcMax", 0.25*H)
gmsh.model.mesh.field.setNumber(6, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(6, "DistMax", r)

gmsh.model.mesh.field.add("Distance",7)
gmsh.model.mesh.field.setNumbers(7, "EdgesList", walls)
gmsh.model.mesh.field.add("Threshold",8)
gmsh.model.mesh.field.setNumber(8, "InField", 7)
gmsh.model.mesh.field.setNumber(8, "LcMin", 1e-6)
gmsh.model.mesh.field.setNumber(8, "LcMax", 4  * H)
gmsh.model.mesh.field.setNumber(8, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(8, "DistMax", 0.25 * H)

min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [2, 4, 6, 8])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)

the above code increase the data in the boundary region but the data in other coordinates decreased as it is shown in the mesh figure:

1111
11111589×315 16 KB

I want a uniform data in the area between boundaries

10

See Gmsh tutorial 10;

2.10 t10: Mesh size fields

at

11

Thanks I fixed the mesh size but now I have these error when I want to solve the problem:

RuntimeError: Newton method failed to converge for non-affine geometry

12

Please provide a minimal reproducible example as that error message with no more context doesn’t help.

13

the complete code is:

############################################
#--------------- Imports ------------------#
############################################

import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import tqdm.autonotebook
import matplotlib as mpl

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               create_vector, create_matrix, set_bc)
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities, refine

from dolfinx import plot

from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)

import dolfinx
import viskex
import time
import pylab as pl
from IPython import display
import matplotlib.tri as tri

mesh_comm = MPI.COMM_WORLD

############################################
#--------------- Mesh creation ------------#
############################################
model_rank = 0
gmsh.initialize()

# L = 2.2
# H = 0.41
# c_x = c_y = 0.2
# r = 0.05
L = 3e-4
H = 5e-5
c_x , c_y = (5e-5, 25e-6)
r = 25e-7
gdim = 2

## Creation of the geometry entities (a cylinder within a rectangle)
if mesh_comm.rank == model_rank:
    rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

## The next step is to substract the obstacle from the channel, such that we do not mesh the interior of the circle
if mesh_comm.rank == model_rank:
    fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
    gmsh.model.occ.synchronize()

## To get GMSH to mesh the fluid, we add a physical volume marker
fluid_marker = 1
if mesh_comm.rank == model_rank:
    volumes = gmsh.model.getEntities(dim=gdim)
    assert (len(volumes) == 1)
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
    gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

## To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side)
## with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by computing the center of mass for each geometrical
## entity.

inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, H / 2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H / 2, 0]):
            outflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
            walls.append(boundary[1])
        else:
            obstacle.append(boundary[1])
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

res_min = 5e-6 #for micro
# res_min =r/1.5
gmsh.model.mesh.field.add("Distance",1)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", obstacle)
gmsh.model.mesh.field.add("Threshold",2)
gmsh.model.mesh.field.setNumber(2, "InField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", 5e-7)
gmsh.model.mesh.field.setNumber(2, "LcMax", 10 * H)
gmsh.model.mesh.field.setNumber(2, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(2, "DistMax",  res_min)

gmsh.model.mesh.field.add("Distance",3)
gmsh.model.mesh.field.setNumbers(3, "EdgesList", outflow)
gmsh.model.mesh.field.add("Threshold",4)
gmsh.model.mesh.field.setNumber(4, "InField", 3)
gmsh.model.mesh.field.setNumber(4, "LcMin", r/3)
gmsh.model.mesh.field.setNumber(4, "LcMax", 6e-6)
gmsh.model.mesh.field.setNumber(4, "DistMin", r)
gmsh.model.mesh.field.setNumber(4, "DistMax", r)

gmsh.model.mesh.field.add("Distance",5)
gmsh.model.mesh.field.setNumbers(5, "EdgesList", inflow)
gmsh.model.mesh.field.add("Threshold",6)
gmsh.model.mesh.field.setNumber(6, "InField", 5)
gmsh.model.mesh.field.setNumber(6, "LcMin", r/3)
gmsh.model.mesh.field.setNumber(6, "LcMax", 0.25*H)
gmsh.model.mesh.field.setNumber(6, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(6, "DistMax", r)

gmsh.model.mesh.field.add("Distance",7)
gmsh.model.mesh.field.setNumbers(7, "EdgesList", walls)
gmsh.model.mesh.field.add("Threshold",8)
gmsh.model.mesh.field.setNumber(8, "InField", 7)
gmsh.model.mesh.field.setNumber(8, "LcMin", 3e-6)
gmsh.model.mesh.field.setNumber(8, "LcMax", 4  * H)
gmsh.model.mesh.field.setNumber(8, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(8, "DistMax", 0.25 * H)

min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [2, 4, 6, 8])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)

### GENERATING THE MESH
## We are now ready to generate the mesh. However, we have to decide if our mesh should
## consist of triangles or quadrilaterals. Here we use second order quadrilateral elements.

if mesh_comm.rank == model_rank:
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 1)
    gmsh.option.setNumber("Mesh.RecombineAll", 1)
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
    gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary",0)
    gmsh.option.setNumber("Mesh.MeshSizeFromPoints",0)
    gmsh.option.setNumber("Mesh.MeshSizeFromCurvature",0)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")

#=============== Loading mesh and boundary markers============================
## As we have generated the mesh, we now need to load the mesh and corresponding
## facet markers into DOLFINx. See https://jsdokken.com/src/tutorial_gmsh.html for explanations
mesh, c, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
mesh.name = "Cylinder"
c.name = "Cell markers"
ft.name = "Facet markers"
# Physical and discretization parameters

t = 0
T = 3                      # Final time
dt =1e-7 #1 / 80                # Time step size
num_steps = 50 #300#int(T / dt)
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(0.001))  # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1000))     # Density
#=================== Boundary conditions===========================

v_cg2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

fdim = mesh.topology.dim - 1

# Define boundary conditions
# U_0 = 1.5 # longitunal/lenghtwise velocity at coordinate (x,y)=(0,H/2)

class InletVelocity():
    def __init__(self, t):
        self.t = t

    def __call__(self, x):
        U_0 = 0.03
        values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = U_0*np.ones((1, x.shape[1]), dtype=PETSc.ScalarType)
        return values
    

# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity(t)
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)))
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V)
# Obstacle
bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V)
bcu = [bcu_inflow, bcu_obstacle, bcu_walls]
# Outlet
bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q)
bcp = [bcp_outlet]

#============== Variational form============================
u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)  #uold
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)

# variational formulation for the first step
f = Constant(mesh, PETSc.ScalarType((0, 0)))
F1 = rho / k * dot(u - u_n, v) * dx
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mu * inner(grad(u + u_n), grad(v)) * dx - dot(p_, div(v)) * dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1)

#  the second step
a2 = form(dot(grad(p), grad(q)) * dx)
L2 = form(-rho / k * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# create the last step
a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u_s, v) * dx - k * dot(nabla_grad(phi), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.JACOBI)

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.MINRES)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

#============ Verification of the implementation compute known physical quantities==============
n = -FacetNormal(mesh)  # Normal pointing out of obstacle
dObs = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle_marker)
u_t = inner(as_vector((n[1], -n[0])), u_)
# drag = form(2 / 0.1 * (mu / rho * inner(grad(u_t), n) * n[1] - p_ * n[0]) * dObs)
# lift = form(-2 / 0.1 * (mu / rho * inner(grad(u_t), n) * n[0] + p_ * n[1]) * dObs)
# micro
drag = form(2 / 2.0e-7 * (mu / rho * inner(grad(u_t), n) * n[1] - p_ * n[0]) * dObs)
lift = form(-2 / 2.0e-7 * (mu / rho * inner(grad(u_t), n) * n[0] + p_ * n[1]) * dObs)
if mesh.comm.rank == 0:
    C_D = np.zeros(num_steps, dtype=PETSc.ScalarType)
    C_L = np.zeros(num_steps, dtype=PETSc.ScalarType)
    t_u = np.zeros(num_steps, dtype=np.float64)
    t_p = np.zeros(num_steps, dtype=np.float64)

tree = bb_tree(mesh, mesh.geometry.dim)
points = np.array([[2.25e-05, 5e-5, 0], [2.75e-05, 5e-5, 0]]) #micro
# points = np.array([[0.15, 0.2, 0], [0.25, 0.2, 0]])
cell_candidates = compute_collisions_points(tree, points)
colliding_cells = compute_colliding_cells(mesh, cell_candidates, points)
front_cells = colliding_cells.links(0)
back_cells = colliding_cells.links(1)
if mesh.comm.rank == 0:
    p_diff = np.zeros(num_steps, dtype=PETSc.ScalarType)

pointsd = mesh.geometry.x

cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidatesd = compute_collisions_points(tree, pointsd)
# Choose one of the cells that contains the point
colliding_cellsd = compute_colliding_cells(mesh, cell_candidatesd, pointsd)
for i, point in enumerate(pointsd):
    if len(colliding_cellsd.links(i)) > 0:
        points_on_proc.append(point)
        cells.append(colliding_cellsd.links(i)[0])

points_on_proc = np.array(points_on_proc, dtype=np.float64)

u_values =np.zeros((num_steps,points_on_proc.shape[0]))
v_values = np.zeros((num_steps,points_on_proc.shape[0]))
p_values = np.zeros((num_steps,points_on_proc.shape[0]))
times = np.zeros(num_steps)

# Solving the time-dependent problem
from pathlib import Path
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, "dfg2D-3-u.bp", [u_], engine="BP4")
vtx_p = VTXWriter(mesh.comm, "dfg2D-3-p.bp", [p_], engine="BP4")
vtx_u.write(t)
vtx_p.write(t)
progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)
for i in range(num_steps):

    progress.update(1)
    # Update current time step
    t += dt
    
    # Update inlet velocity
    inlet_velocity.t = t
    u_inlet.interpolate(inlet_velocity)
    
    # Step 1: Tentative velocity step
    A1.zeroEntries()
    assemble_matrix(A1, a1, bcs=bcu)
    A1.assemble()
    with b1.localForm() as loc:
        loc.set(0)
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_s.vector)
    u_s.x.scatter_forward()

    # Step 2: Pressure corrrection step
    with b2.localForm() as loc:
        loc.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, phi.vector)
    phi.x.scatter_forward()

    p_.vector.axpy(1, phi.vector)
    p_.x.scatter_forward()

    # Step 3: Velocity correction step
    with b3.localForm() as loc:
        loc.set(0)
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)
    u_.x.scatter_forward()

    # Write solutions to file
    vtx_u.write(t)
    vtx_p.write(t)

    times[i]=t
    uval=u_.eval(points_on_proc, cells)
    pval=p_.eval(points_on_proc, cells)
    u_values[i,:]= uval[:,0]
    v_values[i,:]=uval[:,1]
    p_values[i,:] = pval[:,0]

    # Update variable with solution form this time step
    with u_.vector.localForm() as loc_, u_n.vector.localForm() as loc_n, u_n1.vector.localForm() as loc_n1:
        loc_n.copy(loc_n1)
        loc_.copy(loc_n)
        
    
    # Compute physical quantities
    # For this to work in paralell, we gather contributions from all processors
    # to processor zero and sum the contributions.
    drag_coeff = mesh.comm.gather(assemble_scalar(drag), root=0)
    lift_coeff = mesh.comm.gather(assemble_scalar(lift), root=0)
    p_front = None
    if len(front_cells) > 0:
        p_front = p_.eval(points[0], front_cells[:1])
    p_front = mesh.comm.gather(p_front, root=0)
    p_back = None
    if len(back_cells) > 0:
        p_back = p_.eval(points[1], back_cells[:1])
    p_back = mesh.comm.gather(p_back, root=0)
    if mesh.comm.rank == 0:
        t_u[i] = t
        t_p[i] = t - dt / 2
        C_D[i] = sum(drag_coeff)
        C_L[i] = sum(lift_coeff)
        # Choose first pressure that is found from the different processors
        for pressure in p_front:
            if pressure is not None:
                p_diff[i] = pressure[0]
                break
        for pressure in p_back:
            if pressure is not None:
                p_diff[i] -= pressure[0]
                break
vtx_u.close()
vtx_p.close()
t=0

the code works fine with the following mesh:

res_min = 5e-6 
distance_field = gmsh.model.mesh.field.add("Distance")
edges = gmsh.model.getBoundary(volumes, oriented=False)
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", [e[1] for e in edges] )
threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "InField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", r/3)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 10 * H)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", res_min)
min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)
14

The issue is the point evaluation (next time please provide the full stack trace).
I do not have time to work further with this today.
These are the steps that I’ve done to debug the code

for i, (point, cell) in enumerate(zip(points_on_proc, cells)):
        if i == 27:
                uval = u_.eval(point, np.array([cell], dtype=np.int32))

With inspection, I’ve also had a look at the distance from the point to the cell that was determined to be closest:

dolfinx.geometry.compute_distance_gjk(mesh.geometry.x[mesh.geometry.dofmap[cell]], point.reshape(-1, 3))

returning

array([-3.60394061e-07, -6.28160990e-07,  0.00000000e+00])

indicating that the point is ~1e-6 outside the convex hull spanned by the cell in question.

The point in question is:

array([2.99368408e-04, 1.61250761e-06, 0.00000000e+00])

To me it is unclear why you are trying to evaluate the velocity at all nodes of the geometry, as this information is readily available through V.tabulate_dof_coordinates().

15

I tried V.tabulate_dof_coordinates() before but could not update the solution for every timeframes.
So in the current format how I can solve the error?

16

I changed the mesh and the problem is solved.