### Interpolate experimental displacement data on to the mesh surface as dirichlet boundary conditions

106
views
0
10 weeks ago by
Hi friends,

I am solving a hyperelasticty problem to understand the deformation of a material subjected to displacement boundary conditions on the different surfaces of the material. The displacement boundary conditions are obtained from experiments and the format of the data is as follows:

Coordinates:        Displacements
{x,y,z}             {u,v,w}​

The coordinates of the experimental data correspond to arbitrary {x,y,z} coordinates on the surface of the mesh and they do not coincide with nodal coordinates on the surface of the mesh. Is there a way to interpolate this data on to the nodal coordinates of the mesh so that i can apply the data as the displacement boundary conditions? I would appreciate any pointers in this regard.

Thanks

Community: FEniCS Project
You should implement statics steps in your simulation and interpolate previous deformation or transient formulation.
Second you can use Numpy to store your data surface dislocation and interpolate it in your BC (boundary condition) in each step.
written 10 weeks ago by Leo Costa
Thanks @Leo Costa, I don't completely follow. The data that i have is at the final step. Therefore, i will linearly increment the data over a number of static steps and update the the displacement solution from the previous static step. However, i don't quite understand how i would interpolate numpy data on to the boundary and apply it as a boundary condition.
If possible could you please elaborate with sample script on how i would interpolate the stored numpy data and apply it as a boundary condition for a single step. Thanks.
written 10 weeks ago by Kaushik Vijaykumar
Dear All,

I adding the code for the question i mentioned above. I haven't been able figure out how to do it, may be i have not explained the question well. Please find a sample code attached. I would like to apply the boundary condition bct which is tabulated values of displacements for a given x, z coordinates on the boundary of the surface y = 1. The displacement bcs need to be scaled with load steps linearly and the incrementation is performed inside the iteration loop. My apologies for adding such a long code to the community Q&A, as i know its a bad etiquette to do so but i couldn't figure out a way to explain it better. I would really appreciate any help or suggestions in this regard.
Thanks
from dolfin import *
from fenics import *
import numpy as np
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
mesh = UnitCubeMesh(10,10,10)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Boundary identification
#----------------------------------------------------------------------------
# Define boundary sets for boundary conditions
#----------------------------------------------------------------------------
class Left(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-3 # tolerance for coordinate comparisons
return on_boundary and abs(x[0] - 0.0) < tol

class Top(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-3 # tolerance for coordinate comparisons
return on_boundary and abs(x[1] - 1.0) < tol

class Bottom(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-3 # tolerance for coordinate comparisons
return on_boundary and abs(x[1] - 0.0) < tol

class Front(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-3 # tolerance for coordinate comparisons
return on_boundary and abs(x[2] - 0.0) < tol

# Initialize sub-domain instances
left = Left()
top = Top()
bottom = Bottom()
front = Front()

# define meshfunction to identify boundaries by numbers
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1) # mark left as 1
top.mark(boundaries, 2) # mark top as 2
bottom.mark(boundaries, 3) # mark bottom as 3
front.mark(boundaries, 4) # mark front as 4
plot(boundaries,interactive=True)
# Define new measure including boundary naming
ds = Measure("ds")[boundaries] # left: ds(1), right: ds(2), top: ds(3), bottom: ds(4)

bcl = DirichletBC(V.sub(0), Constant(0.0), boundaries, 1)
bcb = DirichletBC(V.sub(1), Constant(0.0), boundaries, 3)
bcf = DirichletBC(V.sub(2), Constant(0.0), boundaries, 4)

# For the top boundary i have a tabluated values of uy displacement for a given x, y and
z coordinates read from SurfaceDisplacementSampleData.npz and stored in data.

# I want to linearly scale the data with 't' for increasing load steps
bct = DirichletBC(V.sub(1), data, boundaries, 2)

bcs = [bcl, bcb, bcf, bct]

load_steps = 3 # number of time steps

# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration

# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
C = F.T*F # Right Cauchy-Green tensor
B = F*F.T

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

# Elasticity parameters
E, nu = 10.0, 0.45
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Stress
stress = (mu/pow(J,5/3))*(B - (1./3.)*tr(B)*I) + (lmbda + (2.*mu/3.))*(J - 1.)*I

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

# Total potential energy
Pi = psi*dx - dot(T, u)*ds(2)

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Non linear variantional problem
problem_u_nl = NonlinearVariationalProblem(F, u, bcs, J)

# Non-linear solvers
solver_u = NonlinearVariationalSolver(problem_u_nl)

# Solver parameters
prm = solver_u.parameters
prm['newton_solver']['linear_solver'] = 'bicgstab'
prm['newton_solver']['preconditioner'] = 'petsc_amg'
prm['newton_solver']['absolute_tolerance'] = 1E-5
prm['newton_solver']['relative_tolerance'] = 1E-5
prm['newton_solver']['maximum_iterations'] = 50
prm['newton_solver']['relaxation_parameter'] = 1.0

# Save solution in VTK format
file = File("u.pvd");
file << u​​