How can I write this equation in a form of EXpression in the fenics


221
views
0
6 months ago by
mamado  
  $T\left(r,t\right)=2\cdot\frac{\left(Ti-T0\right)}{R}\cdot\sum_{i=1}^nJ_0\left(w_n.r\right)\div w_nJ_1\left(w_n.r\right)\cdot e-K\cdot w_nt$T(r,t)=2·(TiT0)R ·i=1nJ0(wn.r)÷wnJ1(wn.r)·eK·wnt


,I tried with this code below  but I don't know how i could do it in the Expression(" "). And I don't know how to use do the summation operator either.

Please help.
from fenics import *
from mshr import *
from scipy.special import j0, j1, jn_zeros
from scipy import linspace
from math import exp
import numpy as np
import numpy as pl
R=1
roots = jn_zeros(0, N) #the first 3 zeros of J0=====> wn
coeff = [1/ (r*j1(r)) for r in roots]
K = 1
X = 2*(Ti-T0)/R
N = 10
for i in range(N) :
F= Expression("X*sum([coeff[i]*j0(roots[i]*(1.-x[0]*x[0]-x[1]*x[1]))*exp(-t*(roots[i]**2))])",t=0.01,degree=1 )
Community: FEniCS Project

1 Answer


1
6 months ago by
Hi mamado,

Try this idea
from fenics import *
from mshr import *
from scipy.special import j0, j1, jn_zeros
from scipy import linspace
from math import exp
import numpy as np
import numpy as pl


mesh = UnitSquareMesh(3, 3)

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

#create a class that inherits form class Expression
class T_func(Expression):
    def __init__(self, R, K, Ti, T0, N,t,degree):
        self.R = R
        self.K = K
        self.Ti = Ti
        self.T0 = T0
        self.N = N
        self.t = t
    def eval(self,value,x):
        Wn = jn_zeros(0, self.N)
        r = 1.0 - x[0]*x[0] - x[1]*x[1] # Is it? or 1 - np.sqrt(x[0]*x[0] + x[1]*x[1])
        sum = 0.0
        for i in range(0,self.N):
            sum = sum + j0(Wn[i]*r)/(Wn[i]*j1(Wn[i]*r))*np.exp(-self.t*Wn[i]**2) 
            #is it? comparing with formulation this line should be 
            #sum = sum + j0(Wn[i]*r)/(Wn[i]*j1(Wn[i]*r))*np.exp(-self.K*Wn[i]*self.t) 

        value[0] = 2*((self.Ti-self.T0)/self.R)*sum

#To use the created class
R = 1
K = 1
N = 10
Ti = 10#your value
T0 = 5#your value
t = 0.01

T_exp = T_func(R=R, K=K, Ti=Ti, T0=T0, N=N,t=t,degree=1)

#To interpolate the values in your mesh

T = interpolate(T_exp, V)​

Good luck
Hi ,Ruben Gonzalez  
thanks foryou. Of course this was quite obvious ...
Best regards,
Mamadou
written 5 months ago by mamado  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »