How to solve an equation with an array data

General
1

Hello everyone,

I am currently attempting to solve an equation using an array of data. Unfortunately, I have encountered an error in the process. I would greatly appreciate your assistance in resolving this matter. Thank you in advance for your help.

from fenics import *
import numpy as np


nx  = 30   # number of element
x_L = 13   # width
a   = x_L / 6
mesh    = IntervalMesh(nx, -x_L, x_L)
domains = MeshFunction("size_t", mesh, 1)
domains.set_all(0)

def boundary(x):
    return x[0] == x_L or x[0] == -x_L
V = FunctionSpace(mesh, "Lagrange", 1)

################################# example of array ##################
F_z = interpolate(Expression(" x[0] ", degree=1), V)
F_zz = np.array(F_z.vector())
data = []
for i in F_zz :
    if i >=  -a/2 and i <= a / 2:
        data.append(0)
    else :
        data.append(1)
data_1 = np.array(data)
##################################### Dirichlet condition ##############
bc = DirichletBC(V, 0.0, boundary)
u = TrialFunction(V)
v = TestFunction(V)
dx = Measure("dx")[domains]
a =  inner(nabla_grad(u), nabla_grad(v)) * dx(1)  +  data_1 * u * v * dx
b =  u * v * dx
# Assemble matrices
H = PETScMatrix()
M = PETScMatrix()
assemble(a, tensor=H)
assemble(b, tensor=M)

######## Apply boundary condition #######################################
bc.apply(H)
bc.apply(M)
# Create eigensolver for generalized EVP
eigensolver = SLEPcEigenSolver(H, M)
eigensolver.parameters["spectrum"] = "smallest magnitude"
eigensolver.parameters["solver"] = "lapack"
# Compute generalized eigenvalue decomposition
print("Computing eigenvalue decomposition. This may take a while.")
eigensolver.solve()
i = 0
E0 = []
while i < eigensolver.get_number_converged():
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    E0.append(r)
    i = i + 1

print("sol = " , E0[0])
2

Hi, I think instead of using numpy array you should use PETSc, therefore please consider reading Conversion from array to PETSc and from PETSc to array.

1 Like
3

Thank you very much for your suggestion. Unfortunately, I tried implementing it, but I didn’t manage to find a solution.

4

Why are you translating a dolfin Function into an array?

Are you trying to make an indicator function for only integrating over a certain subdomain?

1 Like
5

Thank your Mr. dokken for your response. This array is just a sample example. At every position in the system, I have a value that isn’t only related to the subdomain.

6

It has to be related in some way right? How should the vector be multiplied with u and v?

1 Like
7

I solve this matter by using this code:

data = []
for i in F_zz :
    if i >=  -a/2 and i <= a / 2:
        data.append(0)
    else :
        data.append(1)
data_1 = Function(V)
data_1.vector()[:] = data