assemble vector in PointSource function


117
views
0
8 weeks ago by
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}$·σ=i,jPn(x)δ(xxi,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}$xi,j is the midpoint of  $x_i$xi and $x_j$xj,  $\Gamma_{\left\{i,j\right\}}$Γ{i,j} is the distance between  $x_i$xi and  $x_j$xj .

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
8 weeks ago by
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:

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.
Thanks a lot for your kind help! It does work and my plot looks reasonable. Thanks again for helping me out of this suck :)
written 8 weeks ago by Alice Peng  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »