Hello Everyone, using help from the solution of my previous question in the forum Circular Symmetry Boundary Condition on a sector of circle - #7 by dokken.
I tried limiting the rigid body motion occurring in my solution using the rigid motion null-space. However, the displacements are coming random. Expected uy displacement is
Here is the MWE:
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") = {19};
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;
Field[1].YMin = 0;
Field[1].ZMax = 0;
Field[1].ZMin = 0;
Here is the 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
import dolfinx.cpp as _cpp
from dolfinx import mesh, fem, io, common, default_scalar_type, la
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
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]]])
def rigid_motions_nullspace(V: fem.FunctionSpaceBase):
"""
Function to build nullspace for 2D/3D elasticity.
Args:
V: The function space
"""
_x = fem.Function(V)
# Get geometric dim
gdim = V.mesh.geometry.dim
assert gdim == 2 or gdim == 3
# Set dimension of nullspace
dim = 3 if gdim == 2 else 6
# Create list of vectors for null space
nullspace_basis = [_x.vector.copy() for _ in range(dim)]
from contextlib import ExitStack
with ExitStack() as stack:
vec_local = [stack.enter_context(x.localForm()) for x in nullspace_basis]
basis = [np.asarray(x) for x in vec_local]
dofs = [V.sub(i).dofmap.list.reshape(-1) for i in range(gdim)]
# Build translational null space basis
for i in range(gdim):
basis[i][dofs[i]] = 1.0
# Build rotational null space basis
x = V.tabulate_dof_coordinates()
dofs_block = V.dofmap.list.reshape(-1)
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
if gdim == 2:
basis[2][dofs[0]] = -x1
basis[2][dofs[1]] = x0
elif gdim == 3:
basis[3][dofs[0]] = -x1
basis[3][dofs[1]] = x0
basis[4][dofs[0]] = x2
basis[4][dofs[2]] = -x0
basis[5][dofs[2]] = x1
basis[5][dofs[1]] = -x2
la.orthonormalize(nullspace_basis)
assert la.is_orthonormal(nullspace_basis, float(500 * np.finfo(_x.x.array.dtype).resolution))
return PETSc.NullSpace().create(vectors=nullspace_basis) # type: ignore
# 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)
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)
# Dirichlet bc at a point with both ux = 0
def point1(x):
return np.isclose(x[0], 0.05) & np.isclose(x[1], 0)
def point2(x):
return np.isclose(x[0], 0.05) & np.isclose(x[1], 0.1)
constrained_vertex1 = locate_entities(domain, 0, point1)
constrained_vertex2 = locate_entities(domain, 0, point2)
constrained_dof1 = fem.locate_dofs_topological(V.sub(0), 0, constrained_vertex1)
constrained_dof2 = fem.locate_dofs_topological(V.sub(0), 0, constrained_vertex2)
u1_bc = fem.Constant(domain, 0.0)
bc_x1 = fem.dirichletbc(u1_bc, constrained_dof1, V.sub(0))
u2_bc = fem.Constant(domain, 0.0)
bc_x2 = fem.dirichletbc(u2_bc, constrained_dof2, V.sub(0))
bcs = [bc_x1, bc_x2]
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()
print(len(mpc.slaves))
null_space = rigid_motions_nullspace(mpc.function_space)
ksp_rtol = 1e-8 * np.finfo(default_scalar_type).resolution
petsc_options = {
"ksp_rtol": ksp_rtol,
"pc_type": "gamg",
"pc_gamg_type": "agg",
"pc_gamg_coarse_eq_limit": 1000,
"pc_gamg_sym_graph": True,
"mg_levels_ksp_type": "chebyshev",
"mg_levels_pc_type": "jacobi",
"mg_levels_esteig_ksp_type": "cg",
"matptap_via": "scalable",
"pc_gamg_square_graph": 2,
"pc_gamg_threshold": 0.02,
# ,"help": None, "ksp_view": None
}
uh = fem.Function(mpc.function_space, name="Displacement")
WPsolve = LinearProblem(a_form, L_form, mpc, u=uh, bcs=bcs, petsc_options=petsc_options)
WPsolve.A.setNearNullSpace(null_space)
uh = WPsolve.solve()
uh.name = "u"
null_space.remove(uh.vector)
mpc.homogenize(uh)
mpc.backsubstitution(uh)
uh.x.scatter_forward()
