Hi Everyone, I have a question in regards to application of circular periodic boundary condition on the sector of a circle for a plane strain case.
In the following problem definition shown below, I apply a constant force on top edge and periodic boundary condition on side edges. In order to obtain a finite element solution for linear elastic analysis, are described boundary conditions enough or do I need to constrain by applying dirichlet boundary condition some where else also?
As I tried this problem in fenicsx and I am getting large displacements as if the there is no constraint applied. Here is my code.
Point(1) = {0, 0, 0, 0.001};
//+
Point(2) = {0.2, 0.5, 0, 0.001};
//+
Point(3) = {-0.2, 0.5, 0, 0.001};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 1};
//+
Curve Loop(1) = {2, 3, 1};
//+
Plane Surface(1) = {1};
//+
Physical Curve("right", 1) = {1};
//+
Physical Curve("top", 2) = {2};
//+
Physical Curve("left", 3) = {3};
//+
Physical Surface("my_surface", 4) = {1};
mesh_path = "finalb.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(mesh_path, MPI.COMM_WORLD, gdim=2)
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, locate_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)
from dolfinx_mpc.utils import (create_point_to_point_constraint, determine_closest_block, rigid_motions_nullspace)
from dolfinx_mpc import LinearProblem, MultiPointConstraint
gdim = 2
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):
out_x = np.zeros(x.shape)
out_x[0] = x[1]*(-0.2/0.5)
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
# Applying a uniform pressure
T = fem.Constant(domain, 1000000.0)
#Self-weight on the surface
n = FacetNormal(domain)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
L_form = dot(T*n,u_) * ds(2)
bcs = []
with common.Timer("~~Periodic: Compute mpc condition"):
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V.sub(0), facet_markers, 1, 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"
