How to impose a pressure distribution from a data file as a boundary condition (loading)?

1

Hello,

We want to impose a pressure distribution (in space) from a data file as a boundary condition. We already performed simulations with a non-uniform loading (in space or in time), but only with an analytic “Expression”. This “Expression” was then used in the computation of the work of external forces and everything worked well.

For the moment, we can read the data file and save the pressure in a vector, but the problem is that we can’t use it in the computation of the work of external forces as it can not be defined as an “Expression”.

Any help would be appreciated. Thanks in advance!

Lucas

2

are the point data from file aligned with the degrees of freedom?

If so, does the input data have both the coordinate and value associated with it?

If so, you can tabulate the dof coordinates of the relevant function space and match the coordinates with the input and assign it directly to the function vector.

3

Thank you for your quick response,

Yes the point data are aligned with the degrees of freedom and I have chosen an appropriate mesh to avoid interpolating the pressure distribution on the present mesh (for the moment).

Ah yes I think that the part that I am missing is “assign it to the function vector”. My code is inspired from this tutorial:
https://fenicsproject.org/olddocs/dolfin/latest/python/demos/elastodynamics/demo_elastodynamics.py.html

In this tutorial, the loading is uniform in space so it can be defined as an expression:

# Define the loading as an expression depending on t
p = Expression(("0", "t <= tc ? p0*t/tc : 0", "0"), t=0, tc=cutoff_Tc, p0=p0, degree=0)

And then the work of external forces can be computed as:

# Work of external forces
def Wext(u_):
    return dot(u_, p)*dss(3)

In my case, I need that p represents the pressure loading at the boundary, obtained from the data file. So first, I read data and I store it in a list. The problem is that I don’t know how to use this list to compute the work of external forces because in:

# Work of external forces
def Wext(u_):
   return dot(u_, p)*dss(3)

, the dimension of p is equal to the degrees of freedom, I don’t know how to compute it for each facet of the boundary.

Sorry if I am not very clear, I am very new with Fenics and not very used to using FEM.

Lucas

4

Do you have an example or a part of a tutorial that shows how to do this?

Thank you in advance.

5

Consider the following:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(1, 1)

V = FunctionSpace(mesh, "CG", 1)


pressure_data = np.array([0.1, 0.3, 0.4, 0.8])
pressure_coordinates = np.array([[1., 1.], [0, 0], [0., 1.], [1., 0]])
ind2 = np.lexsort((pressure_coordinates[:, 1], pressure_coordinates[:, 0]))


p = Function(V)
dof_coords = V.tabulate_dof_coordinates()
ind = np.lexsort((dof_coords[:, 1], dof_coords[:, 0]))

for (in_index, p_index) in zip(ind2, ind):
    p.vector()[p_index] = pressure_data[in_index]

with XDMFFile("p.xdmf") as xdmf:
    xdmf.write(p)
1 Like
6

Thank you a lot for your help, this solution worked very well !

7

In my case, the pressure distribution is a loading only defined on the left boundary of the rectangle. Instead of working with pressure_data = np.array([0.1, 0.3, 0.4, 0.8]) defined on the whole domain, is it possible to define it only on the left boundary? and use it in :

# Work of external forces
def Wext(u_):
   return dot(u_, p)*dss(3)

as mentionned above

8

You would need to create a mapping between the degrees of freedom on that boundary (their location in the p.vector() and the input data from file.

9

So i adapted what we did before with the degrees of freedom but only for the boundary instead of the whole domain and i get the p vector that need with:

path = 'SCB_outputs/SCB_wall_pressure_double00001.data'
p_SCB = np.loadtxt(path, usecols=1, unpack='True')
y_SCB = np.loadtxt(path, usecols=0, unpack='True')

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
leftboundary = LeftBoundary()

boundarymesh = BoundaryMesh(mesh, 'exterior')
left_mesh = SubMesh(boundarymesh, leftboundary)

W = FunctionSpace(left_mesh, "CG", 1)
p = Function(W, name='loading')

ind2 = [iind2 for iind2 in range (len(y_SCB))]
dof_coords = W.tabulate_dof_coordinates()
ind = np.lexsort((dof_coords[:, 1], dof_coords[:, 0]))

for (in_index, p_index) in zip(ind2, ind):
    print(in_index,p_SCB[in_index])
    p.vector()[p_index] = p_SCB[in_index]

Now the problem is that i am not sure about this :

# Work of external forces
def Wext(u_):
    return u_[0]*p*dss(3)

In fact u_ is defined as u_ = TestFunction(V) with V = VectorFunctionSpace(mesh, "CG", 1)
I have the following errors:

    domains = extract_domains(integrand)
    return sorted(join_domains(domainlist))
TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'

Maybe that it can only work with p being on defined V? because u_ is defined on V?

10

This code is not reproducible, i.e.
I cannot run this code and reproduce the error you are getting.

First of all because you are using a non built in mesh to illustrate the issue.
Secondly, you are not defining all variables in your snippet.

12

Yes sorry here is a part of my code and it reproduces the error that I was mentioning :

from dolfin import *
import numpy as np


# Define mesh
nx = 100
ny = 2500
mesh = RectangleMesh(Point(0., 0.), Point(0.002, 0.001), nx, ny, "crossed") #2D

# Reading input data file composed by y position on the boundary and associated pressure
path = 'SCB_outputs/SCB_wall_pressure_double00001.data'
p_SCB = np.loadtxt(path, usecols=1, unpack='True')
y_SCB = np.loadtxt(path, usecols=0, unpack='True')

def left(x, on_boundary):
    return near(x[0], 0.) and on_boundary

#Subdomain creation to then define the boundary step size dss
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
force_boundary = AutoSubDomain(left)
force_boundary.mark(boundary_subdomains, 3)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

#Necessary to define W functionspace associated to the left boundary
leftboundary = LeftBoundary()
boundarymesh = BoundaryMesh(mesh, 'exterior')
left_mesh = SubMesh(boundarymesh, leftboundary)

V = VectorFunctionSpace(mesh, "CG", 1)
W = FunctionSpace(left_mesh, "CG", 1)

u_ = TestFunction(V)
p = Function(W, name='loading')

# vertical position indices
ind2 = [iind2 for iind2 in range (len(y_SCB))]

dof_coords = W.tabulate_dof_coordinates()
# W degrees of freedom indices
ind = np.lexsort((dof_coords[:, 1], dof_coords[:, 0]))

for (in_index, p_index) in zip(ind2, ind):
    p.vector()[p_index] = p_SCB[in_index]

# Define measure for boundary condition integral
dss = ds(subdomain_data=boundary_subdomains)

# Work of external forces
def Wext(u_):
    return u_[0]*p*dss(3)

error = Wext(u_) #error occurs when computing Wext(u_)

I read a data file to determine the pressure load on the left boundary. This data file is composed of 2500 lines. That is why ny = 2500 in order to first avoid to interpolate the input datas on the mesh that I am using with fenics but you can reproduce the error with any data file and replacing ny by the lines number of the data file.

(I get the following error:

Traceback (most recent call last):
  File "dokken.py", line 56, in <module>
    error = Wext(u_) #error occurs when calculating Wext(u_)
  File "dokken.py", line 54, in Wext
    return u_[0]*p*dss(3)
  File "/home/menezl/anaconda3/envs/fenics/lib/python3.8/site-packages/ufl/measure.py", line 437, in __rmul__
    domains = extract_domains(integrand)
  File "/home/menezl/anaconda3/envs/fenics/lib/python3.8/site-packages/ufl/domain.py", line 355, in extract_domains
    return sorted(join_domains(domainlist))
TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'

)

13

Here you are mixing BoundaryMesh, and SubMesh with your normal mesh. You cannot have a variational form with the combination of the two.
I think you should consider MeshView in @cdaversin framework: GitHub - cdaversin/mixed-dimensional-examples: Code to reproduce numerical examples presented in "Abstractions and automated algorithms for mixed-dimensional finite element methods" (2019)