How to apply Dirichlet BC to a mesh that moves? (and move mesh through applying displacement to every point)


274
views
3
6 months ago by

Hello all,

I'm struggling to correctly apply Dirichlet BC to a mesh that moves. I have a solution on a static reference mesh. I want to use one boundary of this solution (for instance the values on the top of a rectangle) as a Dirichlet BC on a separate mesh whose coordinates change in time.

I have several test cases in the attached code. The first applies the reference mesh solution to the boundary of a separate mesh. This happens successfully if the separate mesh has the same x coordinates in this example, but fails if it is moved in both x and y.

If mesh_1 is moved by editing  the coordinates the DirichletBC fails.

If instead I move the separate mesh with ALE DirichletBC succeeds provided I only move in y.

Ultimately I would like to move a mesh through applying a displacement to every coordinate within the mesh (not just the boundary as with my current ALE implementation) before accurately applying boundary conditions to the moved mesh.

Editing .coordinates() of a mesh achieves the first part of this but the boundary condition is proving problematic.

Any help would be much appreciated!

from fenics import *
import numpy as np

# reference mesh that does not change in time.
mesh_ref = RectangleMesh(Point(0.0,0.0),Point(1.0,1.0),10,10)

# two test meshes, both offset from the initial mesh.
mesh_1 = RectangleMesh(Point(0.0,2.0),Point(1.0,3.0),10,10)
mesh_2 = RectangleMesh(Point(0.0,2.0),Point(1.0,3.0),10,10)

# Set up a vector function space on each mesh
V_ref = VectorFunctionSpace(mesh_ref, "CG", 1)
V1 = VectorFunctionSpace(mesh_1, "CG", 1)
V2 = VectorFunctionSpace(mesh_2, "CG", 1)

# simulate some solution on reference mesh
expr = Expression(("cos(x[0])","exp(x[0])"),degree=3)
u_ref = interpolate(expr,V_ref)

# allow extrapolation on that function
u_ref.set_allow_extrapolation(True)

# Want to match values of solution along the top boundary of the reference mesh to
# lower boundary of the separate mesh.

# itdentify boundary on mesh_1
def boundary2(x, on_boundary):
    return near(x[1],2.0)

bc0 = DirichletBC(V1, u_ref, boundary2)

################################################################################
########  Apply BC to separate mesh ##########
################################################################################

# Set up function on mesh_1 and apply BC.
u2 = Function(V1)
u2.vector().zero()
bc0.apply(u2.vector())

# Check error through identifiying nodes on boundary.

# find dofs
dofs_V1 = V_ref.tabulate_dof_coordinates().reshape((V_ref.dim(),-1))
dofs_V2 = V1.tabulate_dof_coordinates().reshape((V1.dim(),-1))
dofs_V3 = V2.tabulate_dof_coordinates().reshape((V2.dim(),-1))

# find indices corresponding to boundaries 1 and 2.
indices_1 = np.where((dofs_V1[:,1] <= 1.0 + DOLFIN_EPS)&(dofs_V1[:,1] >= 1.0 - DOLFIN_EPS))[0]
indices_2 = np.where((dofs_V2[:,1] <= 2.0 + DOLFIN_EPS)&(dofs_V2[:,1] >= 2.0 - DOLFIN_EPS))[0]
indices_3 = np.where((dofs_V3[:,1] <= 2.0 + DOLFIN_EPS)&(dofs_V3[:,1] >= 2.0 - DOLFIN_EPS))[0]

u1_border = u_ref.vector()[indices_1]
u2_border = u2.vector()[indices_2]

print "2 norm mesh 1 and mesh 2 velocities prior to mesh movement:", np.linalg.norm(u1_border - u2_border)/np.linalg.norm(u1_border)
# machine precision achieved.

################################################################################
########  Apply BC to separate mesh with coordinates moved ##########
################################################################################

# Move mesh_1 vertically by 1
for i_coords in range(len(mesh_1.coordinates())):
	mesh_1.coordinates()[i_coords] += [0.0, 1.0]

u3 = Function(V1)
u3.vector().zero()
# bc1.apply(u3.vector()) # poor results
bc1.apply(u3.vector())

u3_border = u3.vector()[indices_2]

print "2 norm mesh 1 and mesh 2 velocities post manual mesh movement:", np.linalg.norm(u1_border - u3_border)/np.linalg.norm(u1_border)

################################################################################
########  Apply BC to separate mesh with coordinates moved in y using ALE ######
################################################################################
# move mesh with ALE

class Interface_1(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 2.0)

interface_1 = Interface_1()

mf = MeshFunction("size_t", mesh_2, 2)
mf.set_all(0)
interface_1.mark(mf,1)

facets = MeshFunction("size_t", mesh_2,1)
facets.set_all(0)
interface_1.mark(facets,1)

# Extract boundary mesh
bmesh = BoundaryMesh(mesh_2, "exterior", True)
for x in bmesh.coordinates():
    # if interface.inside(x, True):
    x[0] += 0.0
    x[1] += 1.0

# update mesh...
ALE.move(mesh_2, bmesh)

bc3 = DirichletBC(V2, u_ref, facets, 1)

u4 = Function(V2)
u4.vector().zero()
bc3.apply(u4.vector()) # achieves zeros no facets matching..
u4_border = u4.vector()[indices_3]

print "2 norm mesh 1 and mesh 2 velocities post ALE y mesh movement:", np.linalg.norm(u1_border - u4_border)/np.linalg.norm(u1_border)

################################################################################
########  ALE with x moved too ######
################################################################################
# y is already moved from previous step

for x in bmesh.coordinates():
    # if interface.inside(x, True):
    x[0] += 1.0
    x[1] += 0.0

ALE.move(mesh_2, bmesh)

bc3 = DirichletBC(V2, u_ref, facets, 1)

u4 = Function(V2)
u4.vector().zero()
bc3.apply(u4.vector()) # achieves zeros no facets matching..
u4_border = u4.vector()[indices_3]

print "2 norm mesh 1 and mesh 2 velocities post ALE x and y mesh movement:", np.linalg.norm(u1_border - u4_border)/np.linalg.norm(u1_border)
Community: FEniCS Project

1 Answer


3
6 months ago by
Maybe this is "cheating", but one way around the problem is to just build the change-of-variables corresponding to the mesh deformation into the variational form rather than actually moving the mesh.  Then BCs can be applied in a convenient, time-independent reference configuration.  Consider the following simple example:

from dolfin import *

mesh = UnitSquareMesh(10,10)

# (including the change-of-variables in the form can cause the automatically-
# computed quadrature degree to blow up unnecessarily)
dx = dx(metadata={'quadrature_degree': 2})

# space for solution
VuE = FiniteElement("Lagrange",mesh.ufl_cell(),1)
Vu = FunctionSpace(mesh,VuE)

# space for mesh deformation
VFE = VectorElement("Lagrange",mesh.ufl_cell(),1)
VF = FunctionSpace(mesh,VFE)

# set some mapping, F, from the mesh to physical space
# (deforms square into parallelogram in this case)
Fexpr = Expression(("x[1]+x[0]","x[1]"),degree=2)
F = interpolate(Fexpr,VF)

# derivatives needed for change-of-variables
DF = grad(F)
J = det(DF)
def gradx(u):
    # (overkill for scalar-valued u)
    n = rank(u)
    ii = indices(n+2)
    # multivariate chain rule:
    # contract over last index of grad(u) and first index of
    # inv(DF) to change variables in derivative; should work for u of arbitrary
    # rank (scalar, vector, tensor)
    return as_tensor(grad(u)[ii[0:n+1]]*inv(DF)[ii[n],ii[n+1]],
                     ii[0:n]+(ii[n+1],))

# pose Poisson problem in physical space
u = TrialFunction(Vu)
v = TestFunction(Vu)

# use gradients w.r.t. physical space (gradx), and include Jacobian determinant
# (J) in integration measure
a = inner(gradx(u),gradx(v))*J*dx

# NOTE: some care must be taken with more complicated Expressions, as "x[0]"
# and "x[1]" refer to coordinates in the UN-deformed mesh, not physical space
L = inner(Constant(1.0),v)*J*dx

# boundary is defined on the UN-deformed mesh
def boundary(x,on_boundary):
    return on_boundary

u = Function(Vu)
solve(a==L,u,[DirichletBC(Vu,Constant(0.0),boundary),])

# Can plot solution on deformed mesh in Paraview by outputting u and d to
# separate files, loading both of those files, applying the "Calculator"
# filter to compute the identity map for each one, selecting both of the
# Calculator results, applying the "Append Attributes" filter, then applying
# "Warp by Vector" using d, and coloring based on u.  
​


EDIT:  For plotting purposes, it would probably make more sense to define a displacement, say, d, as a Function rather than a mapping, F, then replace DF with

DF = grad(d) + Identity(2)

Hey, that's an interesting solution.

I'd like to check my understanding is okay. I intend to move the mesh on which a fluid problem is solved with the incremental pressure correction scheme. I've attached the variational form below, let me know if more context is necessary.

Using your method I would adjust 'dx' to read 'J*dx' and replace  'nabla_grad' in my variational form with 'nabla_gradx' to account for 'DF' in a similar way to what you have with 'gradx'.

Would it work to replace 'ds' with 'J*ds'?

I take it the computed solution, say of the fluid velocity, would be the velocity as if it were computed on a deformed mesh, but associated with undeformed coordinates so boundary condition application is straightforward?

This would mean I could apply velocities directly to the fluid problem, and similarly use the fluid solution to compute a stress on the undeformed boundary. I would then apply gradx in calculating the stress too.

Appreciate your help.

## Variational setup for Incremental Pressure Correction Scheme

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u0) / k, v)*dx \
   + rho*dot(dot(u0, nabla_grad(u0)), v)*dx \
   + inner(sigma(U, p0), epsilon(v))*dx \
   + dot(p0*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p0), nabla_grad(q))*dx - (1/k)*div(u1)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u1, v)*dx - k*dot(nabla_grad(p1 - p0), v)*dx

# Boundary conditions are applied and problem solved in each time step. 

# Subsequently, fluid stress is calculated to use as a boundary condition on an undeformed structure mesh: 

# deviatoric stress tnesor
tau = mu*(grad(u1) + grad(u1).T)

Dim = mesh.topology().dim()
I = Identity(Dim)

# Sigma
sigma_FSI = project(-p1*I + tau, T_space, solver_type = "mumps", \
    form_compiler_parameters = {"cpp_optimize" : True, "representation" : "uflacs"} )
written 6 months ago by felixnewberry  
2
> Using your method I would adjust 'dx' to read 'J*dx' and replace  'nabla_grad' in my variational form with 'nabla_gradx' to account for 'DF' in a similar way to what you have with 'gradx'.

Yes, just keep in mind that nabla_grad() puts the new index FIRST, in contrast to grad(), which puts the new index last (as one would naturally do with the "comma derivative" in index notation), so the definition of nabla_gradx() would be slightly different.

> Would it work to replace 'ds' with 'J*ds'?

No, surface elements transform via Nanson's formula:

https://en.wikipedia.org/wiki/Finite_strain_theory#Transformation_of_a_surface_and_volume_element

(Mind the difference in notation; F in the linked reference, is like my DF.)

> I take it the computed solution, say of the fluid velocity, would be the velocity as if it were computed on a deformed mesh, but associated with undeformed coordinates so boundary condition application is straightforward?

Yes, this sounds right.

> This would mean I could apply velocities directly to the fluid problem, and similarly use the fluid solution to compute a stress on the undeformed boundary. I would then apply gradx in calculating the stress too.

You can compute the Cauchy stress using gradx, but note that multiplying this by the un-deformed normal will not give the correct traction.  The mapping from un-deformed normals to tractions (per unit un-deformed area) is the 1st Piola--Kirchhoff stress, as defined precisely on Wikipedia:

https://en.wikipedia.org/wiki/Stress_measures
written 6 months ago by David Kamensky  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »