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

221

views

0

$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$

,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 )

`T`(`r`,`t`)=2·(`T``i`−`T`0)`R`·∑_{i=1}^{n}`J`_{0}(`w`_{n}.`r`)÷`w`_{n}`J`_{1}(`w`_{n}.`r`)·`e`−`K`·`w`_{n}`t`,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

Hi mamado,

Try this idea

Good luck

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

Please login to add an answer/comment or follow this question.

thanks foryou. Of course this was quite obvious ...

Best regards,

Mamadou