Hello Everyone, I am solving a plane strain periodic boundary condition problem for linear elastic analysis. In order to constraint the geometry, I am thinking of applying dirichlet on a point on the model.
Here, I want to constraint the point in the present example such it is constrained in x-direction, however, it is free to displace in y-direction. Such that during any applied load in y-direction, it is free to displace in y-direction, however, since x-direction is constrained which is at the x=0.05, my model is constraint enough for solution to run.
My question is how do I apply dirichlet boundary condition on a point in just x-direction.
Here, is my MWE for the problem.
Here is my mesh
cl__1 = 0.001;
Point(1) = {0.01, 0.07000000000000001, 0, cl__1};
Point(2) = {0.03, 0.07000000000000001, 0, cl__1};
Point(3) = {0.03, 0.08, 0, cl__1};
Point(4) = {0.01, 0.08, 0, cl__1};
Point(5) = {0.04, 0.07000000000000001, 0, cl__1};
Point(6) = {0.06, 0.07000000000000001, 0, cl__1};
Point(7) = {0.06, 0.08, 0, cl__1};
Point(8) = {0.04, 0.08, 0, cl__1};
Point(9) = {0.07000000000000001, 0.07000000000000001, 0, cl__1};
Point(10) = {0.09, 0.07000000000000001, 0, cl__1};
Point(11) = {0.09, 0.08, 0, cl__1};
Point(12) = {0.07000000000000001, 0.08, 0, cl__1};
Point(13) = {0, 0.06, 0, cl__1};
Point(14) = {0.1, 0.06, 0, cl__1};
Point(15) = {0.1, 0.09, 0, cl__1};
Point(16) = {0, 0.09, 0, cl__1};
Point(17) = {0, 0, 0, cl__1};
Point(18) = {0.1, 0, 0, cl__1};
Point(19) = {0.15, 0.1, 0, cl__1};
Point(20) = {-0.05, 0.1, 0, cl__1};
Line(1) = {4, 1};
Line(2) = {1, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Line(5) = {8, 5};
Line(6) = {5, 6};
Line(7) = {6, 7};
Line(8) = {7, 8};
Line(9) = {12, 9};
Line(10) = {9, 10};
Line(11) = {10, 11};
Line(12) = {11, 12};
Line(13) = {13, 14};
Line(14) = {14, 15};
Line(15) = {15, 16};
Line(16) = {16, 13};
Line(17) = {17, 18};
Line(18) = {18, 19};
Line(19) = {19, 20};
Line(20) = {20, 17};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Curve Loop(3) = {9, 10, 11, 12};
Plane Surface(3) = {3};
Curve Loop(4) = {1, 2, 3, 4, -8, -7, -6, -5, -12, -11, -10, -9, -15, -14, -13, -16};
Plane Surface(4) = {4};
Curve Loop(5) = {16, 13, 14, 15, -20, -19, -18, -17};
Plane Surface(5) = {5};
Physical Line("bottom") = {17};
Physical Line("right") = {18};
Physical Line("top") = {19};
Physical Line("left") = {20};
Physical Line("EF") = {13};
Physical Line("FG") = {14};
Physical Line("GH") = {15};
Physical Line("HE") = {16};
Physical Line("IJ") = {2};
Physical Line("JK") = {3};
Physical Line("KL") = {4};
Physical Line("LI") = {1};
Physical Line("MN") = {6};
Physical Line("NO") = {7};
Physical Line("OP") = {8};
Physical Line("PM") = {5};
Physical Line("QR") = {10};
Physical Line("RS") = {11};
Physical Line("ST") = {12};
Physical Line("TQ") = {9};
Physical Surface("my_surface") = {5};
Physical Surface("EFGH") = {4};
Physical Surface("IJKL") = {1};
Physical Surface("MNOP") = {2};
Physical Surface("QRST") = {3};
Field[1] = Box;
Field[1].Thickness = 0;
Field[1].VIn = 0;
Field[1].VOut = 0;
Field[1].XMax = 0;
Field[1].XMin = 0;
Field[1].YMax = 0;
Here is my code
import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import sys
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx import mesh, fem, io, common, default_scalar_type
import dolfinx.fem.petsc
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities
from dolfinx.plot import vtk_mesh
from dolfinx.io import XDMFFile
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, MultiPointConstraint
# Mesh
mesh_path = "WPright-medium.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(mesh_path, MPI.COMM_WORLD, gdim=2)
import ufl
gdim = 2
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):
out_x = np.zeros(x.shape)
out_x[0] = x[1]*(-0.05/0.1)
out_x[1] = x[1]
return out_x
#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx
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)
V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)
# bc at point
def pointt(x):
return np.logical_and(np.isclose(x[0], 0.05), np.isclose(x[1], 0))
middof = fem.locate_dofs_geometrical(V, pointt)
u_center = fem.Function(V)
u_center.x.array[:] = 0
bc_center = fem.dirichletbc(u_center, middof)
bcs = [bc_center]
with common.Timer("~~Periodic: Compute mpc condition"):
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V.sub(0), facet_markers, 2, periodic_relation, bcs, default_scalar_type(1))
mpc.finalize()
WPsolve = LinearProblem(a_form, L_form, mpc, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = WPsolve.solve()
uh.name = "u"
