How to import csv file data as an expression

1

Hi everyone,

I am using this code as reference.
https://comet-fenics.readthedocs.io/en/latest/demo/elastodynamics/demo_elastodynamics.py.html

In this code, there is a place where I can define the loading depending on time.
Here, a JIT-compiled expression is used.
p0 = 1.
cutoff_Tc = T/5
# 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)

What I want to do is import a csv data file that depends on time something like this
|Time|Value|
|0.002001|0.000934685|
|0.002668|0.000979293|
|0.003335|0.0010239|
|0.004002|0.00106851|
|0.004669|0.0010239|
|0.005336|0.000934685|

I thought I can use How to Import csv data file to FEniCS? this answer but, I’m not sure if this answers my question.
To use a JIT-compiled expression I need to write it in C++ which I haven’t learned much about.

I am happy about any kind of help. Thanks in advance!

2

Have you had a look at interpolation methods from SciPy?

I think

p.assign(sp.interp(t, x_from_csv, y_from_csv))

should be what you are looking for.

It interpolates the value for p at a given time t in the time loop.

3

Hello,
Thank you very much for your advice.

I think it is np.interp. I tried np.interp(t,t,y) where t is the time data and y is the value data.
I am not very sure about what the p.assign part means. It seems like a C++ function but I cannot include this in my code because it is written in python(fenics2018version). On the other side, if I use p.assign in the Expression I cannot use np.interp because it is a python expression.

4

I don’t understand why you insist on using an Expression to deal with the data from the csv-file.

I have attached an example how I dealt with time dependent table based boundary condition in a transient thermal analysis.

import numpy as np
import scipy as sp

import matplotlib.pylab as plt
import matplotlib.tri as tri
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

from mshr import *
from fenics import *

# Spacial discetisation

domain = Polygon([Point(0,-50), Point(50,-50), Point(50,30), Point(30,50), Point(0,50)])
domain -= Circle(Point(0,0), 30)

mesh = generate_mesh(domain, 40)

# ELement type

V = FunctionSpace(mesh, 'P', 1)

# Materialproperties
rho = 1.0
cp = 1.0
kappa = 5

# Initial values

T_Initial = Constant(0.0)

# Boundary values

Fluid_time = np.array([0.0,10,30,40,100])
Fluid_temp = np.array([0.0,100,100,0,0])

T_Fluid = Constant(100)


T_Extern = 0
HTC = 5

# set initial condition

T_n = interpolate(T_Initial, V)

# boundary condition

def boundary(x, on_boundary):
    tol=0.1
    return on_boundary and near(np.hypot(x[0],x[1]), 30.0, tol)

bc = DirichletBC(V, T_Fluid, boundary)
 
# Time stepping settings

t_ges = 100.0            # final time
num_steps = 500     # number of time steps
dt = t_ges / num_steps # time step size
t = 0.0

# Define variational problem

T, v = TrialFunction(V), TestFunction(V)
f = Constant(0)

F = T*v*dx + dt*kappa*dot(grad(T), grad(v))*dx - (T_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

T = Function(V)

# Time stepping

time = np.array([])
P0 = np.array([])
P1 = np.array([])

# vtkfile = File('try.pvd')

while t <= t_ges:
    
    T_Fluid.assign(sp.interp(t,Fluid_time,Fluid_temp))

    # Update current time
    t += dt

    # Compute solution
    solve(a == L, T, bc)

    time = np.append(time,t)
    P0 = np.append(P0,T([30,0]))
    P1 = np.append(P1,T([30,40]))
#    vtkfile << (T, t)

    # Update previous solution
    T_n.assign(T)
1 Like