Hello everyone,
I have a domain as a square with a circular hole at the middle as in the below picture.
Using BoundaryMesh(mesh, "exterior") will give me both the circular and the square edges.
However, I would like to define the boundary mesh only on the edge of the circular hole.
Is this possible without using the MeshView function from the mixed-dimensional branch?
Thank you very much
Phuris
You can try using SubMesh, e.g.
from dolfin import *
from mshr import *
from vtkplotter.dolfin import plot
class MeshPart(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0.8 and x[0] >= 0.2 and x[1] <= 0.8 and x[1] >= 0.2
meshpart = MeshPart()
domain = Rectangle(Point(0,0), Point(1,1)) - Circle(Point(0.5,0.5),0.25)
mesh = generate_mesh(domain, 20)
sb_mesh = SubMesh(BoundaryMesh(mesh, 'exterior'), meshpart)
plot(sb_mesh)
Thank a lot for your help. This method works.
However, I apologize for not clearly state the purpose that I need to apply the boundary mesh for.
I would like to use the boundary mesh (only the circular edge) as a design variable to calculate the sensitivity. The sensitivity part is calculated using dolfin-adjoint.
By using the submesh, I get an error due to the function transfer_from_boundary(h, meshR) from dolfin-adjoint does not allowed me to use the submesh.
Exception has occurred: AttributeError
'OverloadedMesh' object has no attribute 'entity_map'
File "/workspaces/stokes_nomeshdeform.py", line 58, in <module>
h_V = transfer_from_boundary(h, meshR)
The MWE is provided below (which is adapted from this example)
from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
from mshr import *
set_log_level(LogLevel.ERROR)
# Define classes for edges
class Inflow(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
class Walls(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0) or near(x[1], 1.0)
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return x[0]>0.35 and x[0]<0.65 and x[1]>0.35 and x[1]<0.65
# Create mesh
domain = Rectangle(Point(0,0), Point(1,1)) - Circle(Point(0.5,0.5),0.15)
meshR = generate_mesh(domain, 20)
# Define markers for Sbdomians, Dirichlet and Neuman bcs
markers_boundaries = MeshFunction("size_t", meshR, meshR.geometry().dim()-1, 0)
Inflow().mark(markers_boundaries, 1) # left top edge of fluid mesh
Outflow().mark(markers_boundaries, 2) # middle top edge of fluid mesh
Walls().mark(markers_boundaries, 3) # right top edge of fluid mesh
Obstacle().mark(markers_boundaries, 4) # Left edge of fluid mesh
# Create boundary mesh
b_mesh = SubMesh(BoundaryMesh(meshR, 'exterior'), Obstacle())
S_b = VectorFunctionSpace(b_mesh, "CG", 1)
h = Function(S_b, name="Design")
zero = Constant([0]*meshR.geometric_dimension())
# Create function space to represent h in variational form of s
S = VectorFunctionSpace(meshR, "CG", 1)
h_V = transfer_from_boundary(h, meshR)
h_V.rename("Volume extension of h", "")
# Move the mesh
ALE.move(meshR, h_V)
# Set up the primal problem
# create mixed function space
V2 = VectorElement("CG", meshR.ufl_cell(), 2)
S1 = FiniteElement("CG", meshR.ufl_cell(), 1)
VQ = FunctionSpace(meshR, V2*S1)
# Define test and trial functions and the variational form
(u, p) = TrialFunctions(VQ)
(v, q) = TestFunctions(VQ)
a = inner(grad(u), grad(v))*dx - div(u)*q*dx - div(v)*p*dx
l = inner(zero, v)*dx
# Define boundary conditions
(x,y) = SpatialCoordinate(meshR)
g = Expression(("sin(pi*x[1])","0"),degree=2)
bc_inlet = DirichletBC(VQ.sub(0), g, markers_boundaries, 1)
bc_obstacle = DirichletBC(VQ.sub(0), zero , markers_boundaries, 4)
bc_walls = DirichletBC(VQ.sub(0), zero, markers_boundaries, 3)
bcs = [bc_inlet, bc_obstacle, bc_walls]
# solve for the solution
w = Function(VQ, name="Mixed State Solution")
solve(a==l, w, bcs=bcs)
u, p = w.split()
# Define the objective function
J = assemble(inner(grad(u), grad(u))*dx)
Jhat = ReducedFunctional(J, Control(h))
dJdm = Jhat.derivative()If you are using dolfin-adjoint, there is currently only support for boundary meshes. In the future, i might add submeshes. However, as dolfin-adjoint is not in my main Research area any longer, i cannot give you a date for when this might happen.
With the boundary mesh, as you can see in the Stokes example you have referenced, the mesh is over the whole boundary, while deformation is only hopping over the obstacle. You should be able to do this for any part of the boundary using mesh functions when defining the mesh deformation equations.
Thank you very much for your help.
I would like to ask for a further question.
Is there any example that works with moving the internal edges (e.g. interface between two subdomains)?
Then you cannot use the boundary mesh command.
However, this should be quite straightforward using a volume deformation field.
I am not sure if I understand your reply correctly. Does it mean that instead of defining a boundary mesh, I can directly define a function space and a function object on the domain and directly pass this function object into a mesh deformation function?
Yes. Using a function from the original mesh function space will give you the shape sensitivity for all nodes. You can then set up a deformation scheme with this function, where it is for instance fixed at certain boundaries.
Thank you again for your support. This information is very helpful, I will try it.