I want to define a Dirichlet BC for velocity on multiple surfaces as velocity with some magnitude A in the direction normal to the surface. I am not sure how to do so.
Specifically, I have some mesh defined only on the connected pore space mesh_background_domain,, and a mesh function markers_facetsSolidsIndividual that labels the solid-fluid interface like below,
I want to acheive something similar to this post, but I have no x_c. y_c defined to compute the normal vector numerically. I tried these two versions below, but neither works nor really makes sense. Could someone point me to the right direction? Really appreciate your help. I am not able to upload the specific mesh and mesh function I used, but I am adding the following code that is simlar to how I generated my mesh.
To generate mesh and mesh function
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
# Outer rectangular domain
domain = Rectangle(Point(0, 0), Point(20, 10))
# Define 4 irregular polygons with 6–10 vertices
solid_vertices = [
[Point(4, 2), Point(5.2, 2.3), Point(5.5, 3.1), Point(5.1, 3.8), Point(4.3, 4.2), Point(3.5, 3.5), Point(3.6, 2.6)],
[Point(10, 1), Point(10.8, 1.2), Point(11.6, 2.2), Point(11.4, 3.2), Point(10.7, 3.9), Point(9.8, 3.4), Point(9.2, 2.5), Point(9.5, 1.4)],
[Point(15, 5), Point(15.7, 5.8), Point(16.3, 6.7), Point(15.8, 7.6), Point(15, 8), Point(14.2, 7.5), Point(13.7, 6.8), Point(14.1, 5.7), Point(14.6, 5.2)],
[Point(6, 6), Point(6.8, 6.2), Point(7.4, 6.6), Point(7.7, 7.3), Point(7.3, 8.1), Point(6.5, 8.3), Point(5.8, 7.9), Point(5.6, 7.2)]
]
# Subtract holes from the domain
for vertices in solid_vertices:
poly = Polygon(vertices) # âś… vertices is a list of Points (iterable)
domain -= poly
# Generate mesh
mesh = generate_mesh(domain, 64)
class HoleBoundary(SubDomain):
def __init__(self, points, tol=0.05):
super().__init__()
self.edges = []
self.tol = tol
for i in range(len(points)):
p1 = np.array([points[i].x(), points[i].y()])
p2 = np.array([points[(i+1) % len(points)].x(), points[(i+1) % len(points)].y()])
self.edges.append((p1, p2))
def inside(self, x, on_boundary):
if not on_boundary:
return False
px = np.array([x[0], x[1]])
for p1, p2 in self.edges:
# Compute distance from point to edge
v = p2 - p1
w = px - p1
c1 = np.dot(w, v)
if c1 <= 0:
dist = np.linalg.norm(px - p1)
else:
c2 = np.dot(v, v)
if c2 <= c1:
dist = np.linalg.norm(px - p2)
else:
b = c1 / c2
pb = p1 + b * v
dist = np.linalg.norm(px - pb)
if dist < self.tol:
return True
return False
grain_surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for i, vertices in enumerate(solid_vertices, start=1):
subdomain = HoleBoundary(vertices)
subdomain.mark(grain_surface_markers, i)
source_grains = [3,4]
facet_vals = np.sort(np.unique(grain_surface_markers.array()))
for val in facet_vals:
if val in source_grains:
submesh_coordinates = MeshView.create(grain_surface_markers, val).coordinates()
plt.scatter(submesh_coordinates[:, 0], submesh_coordinates[:, 1],s= 1,label=val)
elif val !=1:
submesh_coordinates = MeshView.create(grain_surface_markers, val).coordinates()
plt.scatter(submesh_coordinates[:, 0], submesh_coordinates[:, 1],s= 1,c='grey')
plt.gca().set_aspect("equal")
plt.show()
Version 1: error “UFLException: Integral of type cell cannot contain a ReferenceNormal.”
V = VectorFunctionSpace(mesh_background_domain, 'P', 2)
n = FacetNormal(mesh_background_domain)
g = Constant(1)
u_normal = project(g * n, V)
u_bc_func = Function(V)
u_bc.interpolate(u_on_facets)
bc = DirichletBC(V, u_bc, markers_facetsSolidsIndividual,[10, 20])
Version 2: error “RuntimeError: Unable to compile C++ code with dijitso”
n = FacetNormal(mesh_background_domain)
source_profile = Expression(('(1+t)*n0', '(1+t)*n1'),g=g,n0=n[0], n1=n[1], t=0,degree=2)
bc = DirichletBC(V, source_profile , markers_facetsSolidsIndividual, [10, 20])