Adding Subboundaries: Warning: Found no facets matching domain for boundary condition

mesh
1

Hi, I am trying to divide my boundaries but when I do that I am having a warning. I shared my code below. Thank you!

image


import gmsh
import meshio
from dolfin import *
from dolfin_adjoint import *


# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True,
               "eliminate_zeros": True,
               "precompute_basis_const": True,
               "precompute_ip_const": True}
a = 1
t = 1
lc = 1
h = 1
x_dir, y_dir = 1, 1
scale = 1
normal = Constant((0, 1, 0))

def createGeometryAndMesh(a, t, lc):
    # Clear all models and create a new one
    gmsh.initialize()
    gmsh.clear()
    gmsh.model.add("t3")

    def geopoints(dx, dy, surf_point):
        P1 = gmsh.model.occ.addPoint(0 + dx, dy, 0, lc)
        P2 = gmsh.model.occ.addPoint(t + dx, dy, 0, lc)
        P3 = gmsh.model.occ.addPoint(dx, h + dy, 0, lc)
        P4 = gmsh.model.occ.addPoint(t + dx, h + dy, 0, lc)

        L1 = gmsh.model.occ.addLine(P1, P2)
        L2 = gmsh.model.occ.addLine(P2, P4)
        L3 = gmsh.model.occ.addLine(P4, P3)
        L4 = gmsh.model.occ.addLine(P3, P1)

        Curve1 = gmsh.model.occ.addCurveLoop([L1, L2, L3, L4])
        Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 11 + surf_point)

        return Surface1


    def rectangularbt(dx, dy, surf_point):
        P50 = gmsh.model.occ.addPoint(-0, scale * 0 + h, 0, lc)
        P60 = gmsh.model.occ.addPoint(scale * -0, scale * 0 + h + a, 0, lc)
        P70 = gmsh.model.occ.addPoint(scale * (0 + dx + t), scale * 0 + h + a, 0, lc)
        P80 = gmsh.model.occ.addPoint(scale * (0 + dx + t), scale * 0 + h, 0, lc)

        L50 = gmsh.model.occ.addLine(P50, P80)
        L60 = gmsh.model.occ.addLine(P80, P70)
        L70 = gmsh.model.occ.addLine(P70, P60)
        L80 = gmsh.model.occ.addLine(P60, P50)
        Curve6 = gmsh.model.occ.addCurveLoop([L50, L60, L70, L80])
        Surface6 = gmsh.model.occ.addPlaneSurface([Curve6], 50 + surf_point)

        return Surface6

    # Entity creation operations
    SurfaceA = geopoints(0, 0, 2000)
    SurfaceRb = rectangularbt(x_dir, y_dir, 6000)
    gmsh.model.occ.synchronize()

    frag, frag_map = gmsh.model.occ.fragment(
        [(2, SurfaceRb)],
        [(2, SurfaceA)],
    )

    print("frag :", frag)
    print("frag_map:", frag_map)
    gmsh.model.occ.synchronize()

    # One group for each of the unit squares
    for i in range(1, len(frag_map)):
        print(" [frag_map[i][0][1]]", [frag_map[i][0][1]])
        pg1 = gmsh.model.addPhysicalGroup(2, [frag_map[i][0][1]])
        gmsh.model.setPhysicalName(2, pg1, f"main_geometry{i - 1}")
    # One group for the top and middle rectangles (frag_map[0] and frag_map[1])
    pg2 = gmsh.model.addPhysicalGroup(2, [frag_map[0][0][1]])
    gmsh.model.setPhysicalName(2, pg2, "filled_geometry")
    # Set mesh size at all points
    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
    gmsh.model.occ.synchronize()
    gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
    gmsh.model.mesh.generate()

    gmsh.write("./Result--" + "geoforstl" + ".geo_unrolled")
    gmsh.write("./Result--" + "geo.msh")
    # colors:
    gmsh.option.setNumber("Geometry.PointNumbers", 1)
    gmsh.option.setColor("Geometry.Points", 255, 165, 0)
    gmsh.option.setColor("General.Text", 255, 255, 255)
    gmsh.option.setColor("Mesh.Points", 255, 0, 0)
    rr, gg, bb, aaa = gmsh.option.getColor("Geometry.Points")
    gmsh.option.setColor("Geometry.Surfaces", rr, gg, bb, aaa)
    gmsh.onelab.set("""[
    {"type":"number",
    "name":"Parameters/Twisting angle",
    "values":[90],
    "min":0,
    "max":120,
    "step":1}]""")

    def checkForEvent():
        action = gmsh.onelab.getString("ONELAB/Action")
        if len(action) and action[0] == "check":
            gmsh.onelab.setString("ONELAB/Action", [""])
            createGeometryAndMesh()
            gmsh.graphics.draw()
        return True

    # Refine the mesh
    gmsh.model.mesh.refine()

    gmsh.write("t3.opt");

    gmsh.finalize()

createGeometryAndMesh(a, t, lc)

# Importing mesh from gmsh and defining surface and boundary markers
msh = meshio.read("./Result--" + "geo.msh")

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]

for cell in msh.cells:
    if cell.type == "tetra":
        tetra_cells = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data

triangle_mesh = meshio.Mesh(points=msh.points,
                            cells=[("triangle", triangle_cells)],
                            cell_data={"name_to_read": [triangle_data]})

meshio.write("./Result--" +  "mesh.xdmf",
             triangle_mesh)
meshio.write("mesh.xdmf", triangle_mesh)


# Create mesh and define function space
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
print(mesh.topology().dim() - 1)

File("./Result--" + "MSH.pvd") << mesh
File("./Result--" + "/MSH.pvd") << mesh
File("./Result--" + "MSH2.pvd") << mf
File("./Result--" "/MSH2" ".pvd") << mf
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)


class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)


b_c = Boundary()
b_c.mark(boundary_markers, 3)
File("./Result--" + "MSH3.pvd") << boundary_markers

# Compiling subdomains
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx_A = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=2)
dx_C = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=1)
dx_main = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=3)
V = VectorFunctionSpace(mesh, "CG", 2)  # Lagrange
W_tensor = TensorFunctionSpace(mesh, "CG", 2)
V_tensor = TensorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "R", 0)
# Define functions
du = TrialFunction(V)  # Incremental displacement
v = TestFunction(V)  # Test function
u = Function(V)  # Displacement from previous iteration
u0 = Function(V)  # Initial displacement

def boundary_bot(x, on_boundary):
    tol = 1E-7
    return on_boundary and near(x[1], 0, tol)


bc_bot = DirichletBC(V, [0, 0, 0], boundary_bot)

def boundary_top1(x, on_boundary):
    tol = 1E-7
    return on_boundary and (near(x[1], (a+h), tol) and x[0] > 1.25)

def boundary_top2(x, on_boundary):
    tol = 1E-7
    return on_boundary and (near(x[1], (a+h), tol) and x[0] < 0.25)



disp1 = Constant([0.0, -0.1, 0.0])
bc_top1 = DirichletBC(V, disp1, boundary_top1)
disp2 = Constant([0.0, -0.2 , 0.0])
bc_top2 = DirichletBC(V, disp2, boundary_top2)
bcs = [bc_bot, bc_top1, bc_top2]


boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(2)
AutoSubDomain(boundary_top1).mark(boundary_subdomains, 1)
AutoSubDomain(boundary_top2).mark(boundary_subdomains, 2)
dss = ds(subdomain_data=boundary_subdomains)


def kinematics(u):
    # Kinematics
    d = u.geometric_dimension()  # Space dimension
    I = Identity(d)  # Identity tensor
    F = variable(I + grad(u))  # Deformation gradient
    C = F.T * F  # Right Cauchy-Green tensor

    # Invariants of deformation tensors
    Ic = tr(C)
    J = det(F)

    return Ic, J, F, C, I

# Elasticity parameters
E, nu = 30.5, 0.0  # main
mu, lmbda = Constant(E / (2 * (1 + nu))), Constant(
    E * nu / ((1 + nu) * (1 - 2 * nu)))

disp1.assign(Constant([0.0, -0.1, 0.0]))
disp2.assign(Constant([0.0, -0.1-0.05, 0.0]))
Ic, J, F, C, I= kinematics(u)

psiA = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2
psi_main = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2
psi = psi_main + psiA

# Total potential energy
Pi_main = psi_main * dx_main
PiA = psiA * dx_A
Pi = Pi_main + PiA

dPi = derivative(Pi, u, v)

solve(dPi == 0, u, bcs,
      form_compiler_parameters=ffc_options)
stressA = diff(psiA, F)
stress_main = diff(psi_main, F)
stress = stress_main + stressA

force = dot(stress, normal)
stored_force= assemble((force[1]) * dss(1))
2

These two conditions do not mark a full facet.
Especially, since your mesh goes from 0 to 2 in x direction, with only 4 elements, what facet has all vertices bigger than 1.25?

Note that to debug such problems, you could vastly reduce your problem to something like:

import gmsh
import meshio
from dolfin import *

a = 1
t = 1
lc = 1
h = 1
x_dir, y_dir = 1, 1
scale = 1


def createGeometryAndMesh(a, t, lc):
    # Clear all models and create a new one
    gmsh.initialize()
    gmsh.clear()
    gmsh.model.add("t3")

    def geopoints(dx, dy, surf_point):
        P1 = gmsh.model.occ.addPoint(0 + dx, dy, 0, lc)
        P2 = gmsh.model.occ.addPoint(t + dx, dy, 0, lc)
        P3 = gmsh.model.occ.addPoint(dx, h + dy, 0, lc)
        P4 = gmsh.model.occ.addPoint(t + dx, h + dy, 0, lc)

        L1 = gmsh.model.occ.addLine(P1, P2)
        L2 = gmsh.model.occ.addLine(P2, P4)
        L3 = gmsh.model.occ.addLine(P4, P3)
        L4 = gmsh.model.occ.addLine(P3, P1)

        Curve1 = gmsh.model.occ.addCurveLoop([L1, L2, L3, L4])
        Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 11 + surf_point)

        return Surface1

    def rectangularbt(dx, dy, surf_point):
        P50 = gmsh.model.occ.addPoint(-0, scale * 0 + h, 0, lc)
        P60 = gmsh.model.occ.addPoint(scale * -0, scale * 0 + h + a, 0, lc)
        P70 = gmsh.model.occ.addPoint(scale * (0 + dx + t), scale * 0 + h + a, 0, lc)
        P80 = gmsh.model.occ.addPoint(scale * (0 + dx + t), scale * 0 + h, 0, lc)

        L50 = gmsh.model.occ.addLine(P50, P80)
        L60 = gmsh.model.occ.addLine(P80, P70)
        L70 = gmsh.model.occ.addLine(P70, P60)
        L80 = gmsh.model.occ.addLine(P60, P50)
        Curve6 = gmsh.model.occ.addCurveLoop([L50, L60, L70, L80])
        Surface6 = gmsh.model.occ.addPlaneSurface([Curve6], 50 + surf_point)

        return Surface6

    # Entity creation operations
    SurfaceA = geopoints(0, 0, 2000)
    SurfaceRb = rectangularbt(x_dir, y_dir, 6000)
    gmsh.model.occ.synchronize()

    frag, frag_map = gmsh.model.occ.fragment(
        [(2, SurfaceRb)],
        [(2, SurfaceA)],
    )

    print("frag :", frag)
    print("frag_map:", frag_map)
    gmsh.model.occ.synchronize()

    # One group for each of the unit squares
    for i in range(1, len(frag_map)):
        print(" [frag_map[i][0][1]]", [frag_map[i][0][1]])
        pg1 = gmsh.model.addPhysicalGroup(2, [frag_map[i][0][1]])
        gmsh.model.setPhysicalName(2, pg1, f"main_geometry{i - 1}")
    # One group for the top and middle rectangles (frag_map[0] and frag_map[1])
    pg2 = gmsh.model.addPhysicalGroup(2, [frag_map[0][0][1]])
    gmsh.model.setPhysicalName(2, pg2, "filled_geometry")
    # Set mesh size at all points
    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
    gmsh.model.occ.synchronize()
    gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
    gmsh.model.mesh.generate()

    gmsh.write("./Result--" + "geoforstl" + ".geo_unrolled")
    gmsh.write("./Result--" + "geo.msh")
    # colors:
    gmsh.option.setNumber("Geometry.PointNumbers", 1)
    gmsh.option.setColor("Geometry.Points", 255, 165, 0)
    gmsh.option.setColor("General.Text", 255, 255, 255)
    gmsh.option.setColor("Mesh.Points", 255, 0, 0)
    rr, gg, bb, aaa = gmsh.option.getColor("Geometry.Points")
    gmsh.option.setColor("Geometry.Surfaces", rr, gg, bb, aaa)
    gmsh.onelab.set("""[
    {"type":"number",
    "name":"Parameters/Twisting angle",
    "values":[90],
    "min":0,
    "max":120,
    "step":1}]""")

    def checkForEvent():
        action = gmsh.onelab.getString("ONELAB/Action")
        if len(action) and action[0] == "check":
            gmsh.onelab.setString("ONELAB/Action", [""])
            createGeometryAndMesh()
            gmsh.graphics.draw()
        return True

    # Refine the mesh
    gmsh.model.mesh.refine()

    gmsh.write("t3.opt")

    gmsh.finalize()


createGeometryAndMesh(a, t, lc)

# Importing mesh from gmsh and defining surface and boundary markers
msh = meshio.read("./Result--" + "geo.msh")

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]

for cell in msh.cells:
    if cell.type == "tetra":
        tetra_cells = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data

triangle_mesh = meshio.Mesh(points=msh.points,
                            cells=[("triangle", triangle_cells)],
                            cell_data={"name_to_read": [triangle_data]})

meshio.write("mesh.xdmf", triangle_mesh)


# Create mesh and define function space
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

# Compiling subdomains
V = VectorFunctionSpace(mesh, "CG", 2)  # Lagrange


def boundary_bot(x, on_boundary):
    tol = 1E-7
    return on_boundary and near(x[1], 0, tol)


bc_bot = DirichletBC(V, [0, 0, 0], boundary_bot)


def boundary_top1(x, on_boundary):
    tol = 1E-6
    print(x[0], x[1], a + h, near(x[1], (a + h), tol))
    return on_boundary and (near(x[1], (a + h), tol) and x[0] > 1.25)


def boundary_top2(x, on_boundary):
    tol = 1E-6
    return on_boundary and (near(x[1], (a + h), tol) and x[0] < 0.25)


disp1 = Constant([0.0, -0.1, 0.0])
bc_top1 = DirichletBC(V, disp1, boundary_top1)
disp2 = Constant([0.0, -0.2, 0.0])
bc_top2 = DirichletBC(V, disp2, boundary_top2)
bcs = [bc_bot, bc_top1, bc_top2]
u = Function(V)
bc_top1.apply(u.vector())
print(u.vector().get_local())

Please keep this in mind for further questions.

1 Like
3

Thank you! @dokken
how can I define it fully?

When I use only boundary_top1 as below, I am having something like this:

image
image1171×619 42.1 KB

def boundary_top1(x, on_boundary):
    tol = 1E-7
    return on_boundary and (near(x[1], (a+h), tol) and x[0] > 1.25)


disp1 = Constant([0.0, -0.1, 0.0])
bc_top1 = DirichletBC(V, disp1, boundary_top1)
bcs = [bc_bot, bc_top1]

I also want to add a load to the other half too.

4

Could you make a sketch of what you want to apply where on your domain?

Its unclear to me what you want to apply where

5

Sure, here is the sketch: @dokken

image
image798×392 6.79 KB

6

@dokken I want to apply the displacements [0,-0.1,0] to gamma1 (disp1), and [0,-0.2,0] (disp2) to gamma2 boundary.