### Problem implementing a point load with Taylor-Hood element

165
views
0
4 months ago by

Hello everybody,

I am currently trying to model the problem of an elastic and incompressible plate submitted to a concentrated force on the middle of the top boundary.

To model the concentrated force, I first created and marked a subdomain for the point of application of the force, then I used the point integral dP to implement it in the variational formulation.
Because it is an incompressible material I use the Taylor-Hood element, but it seems to be incompatible with the use of a point integral.

I get indeed the following error:

*** Error:   Unable to assemble form over vertices.
*** Reason:  Expecting test and trial spaces to only have dofs on vertices for point integrals.
*** Where:   This error was encountered inside Assembler.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset:  4c59bbdb45b95db2f07f4e3fd8985c098615527f

Meaning that I cannot use  $P_2\times P_1$P2×P1 element, but    $P_1\times P_{...}$P1×P... . For example, for an element  $P_1\times P_1$P1×P1 the compilation runs without problem.

Do you have any solution to make the point integral compatible with Taylor-Hood element?
Else, do you have a better approach to adress this problem?

You will find below a minimal code reproducing the error.

Haythem

"""
FEniCS program
The model is used to simulate an elastic material under concentrated load.
"""

from __future__ import division
from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt

# Mechanical variables
E = 1; nu = 0.3
lambda_ = E*nu/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))

# Create mesh and define function space
mesh = UnitSquareMesh(50, 50)

# Define function spaces
V = VectorElement('P', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V*Q)

# Define subdomains for boundary conditions
class LeftRight(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and (near(x[0], 0., tol) or near(x[0], 1., tol))

clamped_boundary = LeftRight()

# Define boundary conditions
bc = DirichletBC(W.sub(0), Constant((0, 0)), clamped_boundary)

# Define subdomains for point loads
def __init__(self, xcoords):
self.xcoords = xcoords
SubDomain.__init__(self)
def inside(self, x, on_boundary):
tol = 1E-2
return near(x[1],  1., tol) and near(x[0],  self.xcoords, tol)

force_markers = VertexFunction("size_t", mesh, 0)
force_markers.set_all(0)

# Define the point integral measure
dP = Measure('dP', domain=mesh, subdomain_data=force_markers)

# Define strain and stress
def epsilon(u):

def epsd(u):
return epsilon(u) - 1./2*nabla_div(u)*Identity(2) # deformations planes

def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
d = u.geometric_dimension()  # space dimension

a = (inner(2*mu*epsd(u), epsd(v)) + p*nabla_div(v) + q*nabla_div(u) - 1./lambda_*p*q)*dx

# Solve variational problem
A = assemble(a)
b = assemble(L)
bc.apply(A, b)

# Compute solution
w = Function(W)
solve(A, w.vector(), b)

# Split the mixed solution using deepcopy
# (needed for further computation on coefficient vector)
(u, p) = w.split(True)

# # Split the mixed solution using a shallow copy
(u, p) = w.split()

# Plot solution
plt.figure()
plot(u, title="velocity")
Community: FEniCS Project