### assemble vector in PointSource function

117

views

0

HI everyone,

I am working on a linear elasticity PDE and the right-hand side includes delta function and a vector. Since the PointSource function in FENICS can only work with constant for 'magnitude' parameter. I was wondering if there is any way to solve this and I have tried several different approximations of Delta function, but the result is not good.

In addition, if I understand correctly, when we apply the PointSource to the RHS, actually we add another term in the RHS rather than modify the existing term.

My equation is

$$ $$ $-\nabla\cdot\sigma=\sum_{i,j}P\textbf{n(x)}\delta\left(x-x_{i,j}\right)\bigtriangleup\Gamma_{i,j}$−∇·

with boundary condition

$\frac{\partial\sigma}{\partial\textbf{n}}+\kappa\textbf{u}=0$∂

Here,

Some part of my code is

I am working on a linear elasticity PDE and the right-hand side includes delta function and a vector. Since the PointSource function in FENICS can only work with constant for 'magnitude' parameter. I was wondering if there is any way to solve this and I have tried several different approximations of Delta function, but the result is not good.

In addition, if I understand correctly, when we apply the PointSource to the RHS, actually we add another term in the RHS rather than modify the existing term.

My equation is

$$ $$ $-\nabla\cdot\sigma=\sum_{i,j}P\textbf{n(x)}\delta\left(x-x_{i,j}\right)\bigtriangleup\Gamma_{i,j}$−∇·

`σ`=∑_{i,j}`P`**n(x)**`δ`(`x`−`x`_{i,j})△Γ_{i,j}with boundary condition

$\frac{\partial\sigma}{\partial\textbf{n}}+\kappa\textbf{u}=0$∂

`σ`∂**n**+`κ`**u**=0 .Here,

**n(x)**is unit inward direction and $x_{i,j}$`x`_{i,j}is the midpoint of $x_i$`x`_{i}and $x_j$`x`_{j}, $\Gamma_{\left\{i,j\right\}}$Γ_{{i,j}}is the distance between $x_i$`x`_{i}and $x_j$`x`_{j}.Some part of my code is

```
from __future__ import print_function
from fenics import *
from mshr import *
from dolfin import *
import numpy as np
### define the Euler's Distance of the vectors
def dist(vec1, vec2):
dist = np.linalg.norm(vec1 - vec2)
return dist
# calculate the unit inward direction
def inwardv(center,point1,point2): # center can be circle center/mesh(triangle center)
c=np.array([line_a(point1,point2),line_b(point1,point2)])
midpoint=(point1+point2)/2
if np.dot(c,center-midpoint)>=0:# and np.dot(c,center-midpoint)/(np.linalg.norm(c)*np.linalg.norm(center-midpoint))>=0:
return c/np.linalg.norm(c)
else:
return -c/np.linalg.norm(c)
x0=20
y0=15
nx=40
ny=20
mesh=RectangleMesh(Point(-x0,-y0),Point(x0,y0),nx,ny)
subdomain_vert=np.array([[0,0,1],[0,1.5,1.5]])
center=np.mean(subdomain_vert,axis=1)
### FEM to solve the PDE
### parameters value
E=1
nu=0.49
P=2.08
kappa=3
V = VectorFunctionSpace(mesh, 'P', 1)
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
return E/(1+nu)*(epsilon(u)+tr(epsilon(u))*Identity(d)*nu/(1-2*nu))
u = Function(V)
v = TestFunction(V)
d = u.geometric_dimension() # space dimension
W= TensorFunctionSpace(mesh,'CG',1)
a = inner(sigma(u), grad(v))*dx+kappa*inner(u,v)*ds
L = inner(Constant((0,0)),v)*dx
A,b=assemble_system(a,L)
#### To build up the point source #####
for i in np.arange(subdomain_vert.shape[1]):
point1=subdomain_vert[:,i]
for j in np.arange(i+1,subdomain_vert.shape[1]):
point2=subdomain_vert[:,j]
inwardn = inwardv(center, point1, point2)
midpoint = (point1 + point2) / 2
# b0=assemble(inner(inwardn,v)*dx)
# delta=PointSource(V,Point(midpoint[0], midpoint[1]), P*dist(point1, point2))
# delta.apply(b0)
# b=b+b0
solve(A,u.vector(),b)
```

I has been stuck in this for more than one month, any help will be appreciated.

Best,

Alice

P.S. I have also tried to split the tensor to piecewise the equation, it seems not work either, see https://www.allanswered.com/post/baznk/approach-and-split-the-tensor-component/

Community: FEniCS Project

### 1 Answer

3

I'm having some trouble imagining a problem where you'd want to sum over line segments between every pair of points in a set, but, if you are summing over a collection of line segments in some polyline, it starts to look like numerical quadrature of a traction over an immersed boundary, which is something I've thought about a bit. The following example seems to get the qualitatively-right result and might be helpful:

Note in particular that scalar components of vector point sources can be applied to sub-spaces of the VectorFunctionSpace.
Thanks a lot for your kind help! It does work and my plot looks reasonable. Thanks again for helping me out of this suck :)

```
from dolfin import *
from numpy import array, zeros
import math
from numpy.linalg import norm as npNorm
# Space dimension
d = 2
# Generate a set of n integration elements for a circle
def getCircleSegments(center,r,n):
a = zeros((n,2,d))
dtheta = 2.0*math.pi/float(n)
for i in range(0,n):
theta = i*dtheta
x0 = center[0] + r*math.cos(theta)
y0 = center[1] + r*math.sin(theta)
x1 = center[0] + r*math.cos(theta+dtheta)
y1 = center[1] + r*math.sin(theta+dtheta)
a[i,0,0] = x0
a[i,0,1] = y0
a[i,1,0] = x1
a[i,1,1] = y1
return a
# Generate a flat line at y=0.3 for testing
def getLineSegments(n):
a = zeros((n,2,d))
dx = -1.0/float(n) # go from right to left for downward normal
for i in range(0,n):
x = 1.0+i*dx
x0 = x
y0 = 0.3
x1 = x+dx
y1 = 0.3
a[i,0,0] = x0
a[i,0,1] = y0
a[i,1,0] = x1
a[i,1,1] = y1
return a
# Get the normal to a segment, rotating counter-clockwise by convention
# (will be inward for CCW circle defined above)
def getNormal(p0,p1):
dx = p1[0]-p0[0]
dy = p1[1]-p0[1]
return array([-dy,dx])/npNorm(p1-p0)
# Set up linearized elasticity PDE
mesh = UnitSquareMesh(100,100)
V = VectorFunctionSpace(mesh,"Lagrange",2)
u = TestFunction(V)
v = TrialFunction(V)
def eps(u):
return 0.5*(grad(u) + grad(u).T)
BULK_MOD = Constant(1.0)
SHEAR_MOD = Constant(1.0)
I = Identity(d)
def sigma(u):
return BULK_MOD*tr(eps(u))*I + 2.0*SHEAR_MOD*(eps(u) - tr(eps(u))*I/3.0)
f = Constant((0.0,-1.0))
# Constrain normal displacement on three sides:
def boundary1(x, on_boundary):
return on_boundary and near(x[1],0.0)
bc1 = DirichletBC(V.sub(1), Constant(0.0), boundary1)
def boundary2(x, on_boundary):
return on_boundary and near(x[0],0.0)
bc2 = DirichletBC(V.sub(0), Constant(0.0), boundary2)
def boundary3(x, on_boundary):
return on_boundary and near(x[0],1.0)
bc3 = DirichletBC(V.sub(0), Constant(0.0), boundary3)
# Assemble PDE part
A = assemble(inner(sigma(u),eps(v))*dx)
b = assemble(inner(f,v)*dx)
bc1.apply(A,b)
bc2.apply(A,b)
bc3.apply(A,b)
####### Immersed surface assembly #######
# Magintude of force per unit area applied on the immersed surface
P = 1.0
# Define a series of surface elements
center = array([0.5,0.5])
r = 0.25
n = 1000
segments = getCircleSegments(center,r,n)
#segments = getLineSegments(n)
# For each surface integration element, put a point source in the center,
# with magnitude P*(area of segment) and direction normal to the segment
pointSourceData = [[],[]]
for segment in segments:
p0 = segment[0]
p1 = segment[1]
p = 0.5*(p0+p1)
if(abs(p[0]-0.5)>0.5-DOLFIN_EPS or abs(p[1]-0.5)>0.5-DOLFIN_EPS):
continue
n = getNormal(p0,p1)
w = npNorm(p1-p0) # integration weight
dF = P*w*n
for i in range(0,d):
pointSourceData[i] += [(Point(p[0],p[1]),dF[i]),]
# Apply scalar point sources to sub-spaces corresponding to different
# displacement components
multiPointSources = []
for i in range(0,d):
multiPointSources += [PointSource(V.sub(i),pointSourceData[i]),]
# Apply point sources in batches for each dimension
for i in range(0,d):
multiPointSources[i].apply(b)
#########################################
u = Function(V)
solve(A,u.vector(),b)
# Plot the 11 stress component in the domain; for parts of $\Gamma$ that are
# horizontal, the jump in this should equal $P$ in magnitude
import matplotlib.pyplot as plt
stress = sigma(u)[1,1]
plot(stress)
plt.show()
# PV output
File("stress.pvd") << project(stress,FunctionSpace(mesh,"DG",0))
```

Note in particular that scalar components of vector point sources can be applied to sub-spaces of the VectorFunctionSpace.

written
8 weeks ago by
Alice Peng

Please login to add an answer/comment or follow this question.