Periodic boundary conditions implemented on non-straight Facets in a geometry

dolfinx
, ,
1

In the following geometry, I wish to define periodic boundary conditions on the left and right boundary for a linear elasticity analysis. However, I have hard time implementing it. I have defined the left and right boundary with facet tags. While implementing I get error.

gmsh.initialize()

gmsh.model.add("Winding Pack-1")

gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
h = 1

lc = 1e-3

# Added Points in the geometry

gmsh.model.geo.addPoint(0,0,0,lc,1)
gmsh.model.geo.addPoint(1,0,0,lc,2)
gmsh.model.geo.addPoint(1.5,1,0,lc,3)
gmsh.model.geo.addPoint(-0.5,1,0,lc,4)
#gmsh.model.geo.addPoint(0.3,0.3,0,lc,5)
#gmsh.model.geo.addPoint(0.7,0.3,0,lc,6)
#gmsh.model.geo.addPoint(0.7,0.6,0,lc,7)
#gmsh.model.geo.addPoint(0.3,0.6,0,lc,8)

# Added Lines in the geometry

gmsh.model.geo.addLine(1,2,1)
gmsh.model.geo.addLine(2,3,2)
gmsh.model.geo.addLine(3,4,3)
gmsh.model.geo.addLine(4,1,4)
#gmsh.model.geo.addLine(5,6,5)
#gmsh.model.geo.addLine(6,7,6)
#gmsh.model.geo.addLine(7,8,7)
#gmsh.model.geo.addLine(8,5,8)

gmsh.model.geo.addCurveLoop([1,2,3,4],1)

gmsh.model.geo.addPlaneSurface([1],1)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(1, [1], name="bottom")
gmsh.model.addPhysicalGroup(1, [2], name="right")
gmsh.model.addPhysicalGroup(1, [3], name="top")
gmsh.model.addPhysicalGroup(1, [4], name="left")

gmsh.model.addPhysicalGroup(2, [1], name="my_surface")

gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)

gmsh.model.mesh.generate(gdim)

domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
#facet_markers.name = "Facet markers"

gmsh.write("Winding Pack-1.msh")
gmsh.fltk.run()
gmsh.finalize()

import ufl

def strain(u, repr ="vectorial"):
    eps_t = sym(grad(u))
    if repr =="vectorial":
        return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
    elif repr =="tensorial":
        return eps_t

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
    sigv = dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))


#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform pressure
#rho = 7850
#g = 9.81
T = fem.Constant(domain, 1000000.0)
#fg = fem.Constant(domain, np.array([0, -rho*g]))
#fg = fem.Constant(domain, 78500.0)

#Self-weight on the surface
n = FacetNormal(domain)
#L_form = dot(-fp*n,u_) * ds(5) + dot(fg*n,u_) * ds(5)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
L_form = dot(T*n,u_) * ds(3)
#L_form = dot(-T*n,u_) * ds(5) + dot(fg*n,u_) * ds(5)

#Boundary Conditions (simply supported on top and bottom edges)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))


ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0)),fem.dirichletbc(uy0, top_dofs, V.sub(1)), fem.dirichletbc(ux0, right_dofs, V.sub(0))]
bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0))]



#------------------------------------------------------------------------------------------------------------------------------------------
periodic_boundary = facet_markers.find(2)

#Cyclic Symmetry Condition
def periodic_relation(x):
    out_x = np.zeros(x.shape)
    out_x[0] = 1 - x[0]
    out_x[1] = x[1]
    out_x[2] = x[2]
    return out_x

# Periodic boundary - topological
facets = mesh.locate_entities_boundary(domain, gdim - 1, periodic_boundary)
arg_sort = np.argsort(facet_markers)
mt = mesh.meshtags(domain, gdim - 1, facet_markers[arg_sort], np.full(len(facet_markers), 2, dtype=np.int32))

with common.Timer("~~Periodic: Compute mpc condition"):
    mpc = dolfinx_mpc.MultiPointConstraint(V)
    mpc.create_periodic_constraint_topological(V, mt, 2, periodic_relation, bcs)
    mpc.finalize()

#-------------------------------------------------------------------------------------------------------------------------------------------                      

WPsolve = LinearProblem(a_form, L_form, mpc, u=u, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
WPsolve.solve()

#WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
#WPsolve.solve()


V0 = fem.functionspace(domain, ("DG", 0, (3,)))
sig_exp = fem.Expression(stress(u), V0.element.interpolation_points())
sig = fem.Function(V0, name="Stress")
sig.interpolate(sig_exp)

import pyvista
from dolfinx import plot

pyvista.set_jupyter_backend("static")

topology, cell_types, geometry = plot.vtk_mesh(domain, 2)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

p = pyvista.Plotter()
p.add_mesh(grid,show_edges=True)
p.view_xy()
p.show_axes()
p.show()

vtk = io.VTKFile(domain.comm, "wpex3.pvd", "u")
vtk.write_function(u, 0)
vtk.write_function(sig, 0)
vtk.close()

Problem-statement
Problem-statement859×537 136 KB

2

You have already asked this question at: Implementation of periodic boundary condition for Plane Strain Case- mesh imported using GMSH API
Please remove one of the questions.

3

I have flagged for the removal of this question asked previously. How do I implement the periodic bc for facets 2 (right) and 4 (left) eges as per the geometry along the x-direction.

4

You know the geometric relation for the line on the left and the line on the right.
Thus given appropriate boundary tags of either set of facets, you can make a periodic map

def map(x):
     return (your_map(x[0], x[1]), x[1])

which can be used in the same way as in:

However, as you know that the points are matching, you could also use create_point_to_point_constraint as shown in: dolfinx_mpc/python/demos/demo_elasticity_disconnect.py at 02338e6dcf44f7b5b08238586b277aac8153d4f3 · jorgensd/dolfinx_mpc · GitHub

5

Hello Jorgen,

I did not understand exactly how to create mapping using facets, however, I tried this another approach taking some assistance from your suggested link. I am still getting this error which I am not able to understand. Here is the MWE.

import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import sys

from mpi4py import MPI
from petsc4py import PETSc
#import pyvista

from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
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 meshtags_from_entities, locate_entities_boundary, meshtags
from dolfinx import default_scalar_type

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

gmsh.initialize()

gmsh.model.add("mesh")

gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
h = 0.05

lc = 1e-3

# Added Points in the geometry

gmsh.model.geo.addPoint(0,0,0,lc,1)
gmsh.model.geo.addPoint(0.1,0,0,lc,2)
gmsh.model.geo.addPoint(0.157735027,0.1,0,lc,3)
gmsh.model.geo.addPoint(-0.057735027,0.1,0,lc,4)

# Added Lines in the geometry

gmsh.model.geo.addLine(1,2,1)
gmsh.model.geo.addLine(2,3,2)
gmsh.model.geo.addLine(3,4,3)
gmsh.model.geo.addLine(4,1,4)

gmsh.model.geo.addCurveLoop([1,2,3,4],1)

gmsh.model.geo.addPlaneSurface([1],1)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(1, [1], name="bottom")
gmsh.model.addPhysicalGroup(1, [2], name="right")
gmsh.model.addPhysicalGroup(1, [3], name="top")
gmsh.model.addPhysicalGroup(1, [4], name="left")

gmsh.model.addPhysicalGroup(2, [1], name="my_surface")

gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)

gmsh.model.mesh.generate(gdim)

domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

gmsh.write("mesh.msh")
gmsh.fltk.run()
gmsh.finalize()

import ufl

def strain(u, repr ="vectorial"):
    eps_t = sym(grad(u))
    if repr =="vectorial":
        return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
    elif repr =="tensorial":
        return eps_t

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
    sigv = dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))


#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform traction

T = fem.Constant(domain, 1000000.0)


n = FacetNormal(domain)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
L_form = dot(T*n,u_) * ds(3)

#Boundary Conditions (simply supported on top and bottom edges)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))


ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)

bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]


periodic_boundary = right_dofs  # This is my slave side 

# Periodic boundary - topological
facets = locate_entities_boundary(domain, domain.topology.dim-1, periodic_boundary)
arg_sort = np.argsort(facets)
mt = meshtags(domain, domain.topology.dim -1, facets[arg_sort], np.full(len(facets), 2, dtype=np.int32))

#Cyclic Symmetry Condition
def periodic_relation(x):
    out_x = np.zeros(x.shape)
    out_x[0] = ((np.pi)/3) - x[0]
    out_x[1] = x[1]
    #out_x[2] = x[2]
    return out_x

with common.Timer("~~Periodic: Compute mpc condition"):
    mpc = dolfinx_mpc.MultiPointConstraint(V)
    mpc.create_periodic_constraint_topological(V.sub(0), mt, 2, periodic_relation, bcs, default_scalar_type(1))
    mpc.finalize()
                     

WPsolve = LinearProblem(a_form, L_form, mpc, u=u, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
WPsolve.solve()

My error is

TypeError                                 Traceback (most recent call last)
Cell In[16], line 72
     69 periodic_boundary = right_dofs  # This is my slave side 
     71 # Periodic boundary - topological
---> 72 facets = locate_entities_boundary(domain, domain.topology.dim-1, periodic_boundary)
     73 arg_sort = np.argsort(facets)
     74 mt = meshtags(domain, domain.topology.dim -1, facets[arg_sort], np.full(len(facets), 2, dtype=np.int32))

File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/mesh.py:233, in locate_entities_boundary(mesh, dim, marker)
    208 def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:
    209     """Compute mesh entities that are connected to an owned boundary
    210     facet and satisfy a geometric marking function.
    211 
   (...)
    231 
    232     """
--> 233     return _cpp.mesh.locate_entities_boundary(mesh._cpp_object, dim, marker)

TypeError: locate_entities_boundary(): incompatible function arguments. The following argument types are supported:
    1. (mesh: dolfinx.cpp.mesh.Mesh_float32, dim: int, marker: Callable[[numpy.ndarray[numpy.float32]], numpy.ndarray[bool]]) -> numpy.ndarray[numpy.int32]
    2. (mesh: dolfinx.cpp.mesh.Mesh_float64, dim: int, marker: Callable[[numpy.ndarray[numpy.float64]], numpy.ndarray[bool]]) -> numpy.ndarray[numpy.int32]

Invoked with: <dolfinx.cpp.mesh.Mesh_float64 object at 0x71c5bedc4fb0>, 1, [array([6212, 6218, 6224, 6230, 6236, 6242, 6248, 6262, 6268, 6282, 6288,
       6294, 6300, 6306, 6312, 6318, 6324, 6330, 6336, 6342, 6348, 6354,
       6358, 6360, 6362, 6364, 6366, 6368, 6370, 6372, 6376, 6382, 6384,
       6388, 6392, 6398, 6400, 6402, 6404, 6406, 6408, 6410, 6412, 6414,
       6416, 6418, 6420, 6422, 6424], dtype=int32), array([3106, 3109, 3112, 3115, 3118, 3121, 3124, 3131, 3134, 3141, 3144,
       3147, 3150, 3153, 3156, 3159, 3162, 3165, 3168, 3171, 3174, 3177,
       3179, 3180, 3181, 3182, 3183, 3184, 3185, 3186, 3188, 3191, 3192,
       3194, 3196, 3199, 3200, 3201, 3202, 3203, 3204, 3205, 3206, 3207,
       3208, 3209, 3210, 3211, 3212], dtype=int32)]
6

Here you are mixing dofs and facets:

You have stated that

periodic_boundary = right_dofs  # This is my slave side 

# Periodic boundary - topological
facets = locate_entities_boundary(domain, domain.topology.dim-1, periodic_boundary)
arg_sort = np.argsort(facets)

First of all, periodic_boundary is a list of degrees of freedom.
You are trying to send these into locate_entities_boundary, which takes in a function for locating entities (facets, edges or vertices) based on a geometrical function.

Here is a minimal working example for your case. Please inspect it carefully as I have taken quite a lot of care reducing it to a minimal example

import gmsh
import numpy as np
from mpi4py import MPI

from dolfinx import mesh, fem
from dolfinx.common import Timer
import dolfinx.fem.petsc

from dolfinx.io import VTXWriter, gmshio
from dolfinx import default_scalar_type

from ufl import (FacetNormal, as_matrix, TestFunction, TrialFunction, as_vector, dot, ds, dx, inner, grad, sym)
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem

gmsh.initialize()

gmsh.model.add("mesh")

gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
h = 0.05

lc = 1e-3

# Added Points in the geometry

delta_x = 0.57735027
gmsh.model.geo.addPoint(0,0,0,lc,1)
gmsh.model.geo.addPoint(0.1,0,0,lc,2)
gmsh.model.geo.addPoint(0.1+delta_x,0.1,0,lc,3)
gmsh.model.geo.addPoint(-delta_x,0.1,0,lc,4)

# Added Lines in the geometry

gmsh.model.geo.addLine(1,2,1)
gmsh.model.geo.addLine(2,3,2)
gmsh.model.geo.addLine(3,4,3)
gmsh.model.geo.addLine(4,1,4)

gmsh.model.geo.addCurveLoop([1,2,3,4],1)

gmsh.model.geo.addPlaneSurface([1],1)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(1, [1], name="bottom")
gmsh.model.addPhysicalGroup(1, [2], name="right")
gmsh.model.addPhysicalGroup(1, [3], name="top")
gmsh.model.addPhysicalGroup(1, [4], name="left")

gmsh.model.addPhysicalGroup(2, [1], name="my_surface")

gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)

gmsh.model.mesh.generate(gdim)

domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

gmsh.write("mesh.msh")
#gmsh.fltk.run()
gmsh.finalize()

import ufl

def strain(u, repr ="vectorial"):
    eps_t = sym(grad(u))
    if repr =="vectorial":
        return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
    elif repr =="tensorial":
        return eps_t

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
    sigv = dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))

#Cyclic Symmetry Condition
def periodic_relation(x):
    """
    Map the line 
    y = 0.1/delta_x * x - 0.01/delta_x
    to
    y = -0.1/delta_x*x
    to
    meaning that if we have (x,y) from the right line 
    x_left = -delta_x/0.1 * y
    if y = 0 -> x_left = 0
    if y = 0.01 -> x_left = -delta_x
    """
    out_x = np.zeros(x.shape)
   
    
    out_x[0] = x[1] * (-delta_x/0.1)
    out_x[1] = x[1]
    return out_x

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))

ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)

bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]

with Timer("~~Periodic: Compute mpc condition"):
    mpc = dolfinx_mpc.MultiPointConstraint(V)
    mpc.create_periodic_constraint_topological(V.sub(0), facet_markers, 2, periodic_relation, bcs, default_scalar_type(1))
    mpc.finalize()

#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(mpc.function_space, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform traction

T = fem.Constant(domain, 1000000.0)


n = FacetNormal(domain)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
L_form = dot(T*n,u_) * ds(3)



WPsolve = LinearProblem(a_form, L_form, mpc, u=u, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
WPsolve.solve() 

with dolfinx.io.VTXWriter(domain.comm, "solution.bp", [u]) as bp:
    bp.write(0)
1 Like
7

Many Thanks, Jorgen! It is working now. Now I understand my mistake.