Hello everyone,
I really need your help. I’ve managed to solve a problem with FEniCSx, but I’m having trouble displaying the solution paraview. It seems that the solution is correctly present, as you can see in the attached image.
However, to my surprise, the only result displayed is the one corresponding to the last solution written in the code (ui
). Below, you’ll find a part of the code where I split the solution Ksol
and save it in xdmf format.
I don’t understand why the solution is marked as (partial). I’m not sure if I’m saving the solutions incorrectly.
Thank you in advance for your assistance.
# import pyhton and Fenicsx dependecy
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import geometry,fem, io
from dolfinx.io import gmshio, XDMFFile
from dolfinx.fem import FunctionSpace, Function, Constant, Expression
from dolfinx.fem.petsc import LinearProblem, assemble, NonlinearProblem, assemble_matrix, create_vector, assemble_vector
from dolfinx.nls.petsc import NewtonSolver
import dolfinx.fem.petsc
import ufl
from ufl import dx, grad, inner, ds, Measure, TrialFunctions, TestFunctions, SpatialCoordinate, lhs, rhs, split, lt, gt
........
#import the gmsh generated by gmsh: it should be in gmsh format
msh, cell, facet_tags = gmshio.read_from_msh("GI.msh",MPI.COMM_WORLD, gdim=3)
# time intergration
t = 0.0; dt = 0.1; Tfinal=100
freqSave = 10; inc = 0
stim_t0 = dt; stim_t2 = 5000; stim_dur = 10; stim_amp = 1
# ******* Create mesh and define function spaces **********
P1 = ufl.FiniteElement("CG", msh.ufl_cell(), 1)
P2v = ufl.VectorElement("CG", msh.ufl_cell(), 2)
FiberSpace = FunctionSpace(msh,("CG",1))
Mh = FunctionSpace(msh, P1)
Nh = FunctionSpace(msh, ufl.MixedElement([P1,P1,P1,P1]))
phih = Function(Mh)
#print(" **************** Electr Dof = ", Nh.dim())
# Define Test and Trial function
(ul,vl,ui,vi) = TrialFunctions(Nh)
(w1,w2,m1,m2) = TestFunctions(Nh)
Ksol = Function(Nh)
#ul,vl,ui,vi = ufl.split(Ksol)
#Ksol.name = "Electrophysiology"
Ksol0 = Function(Nh)
ul_old,vl_old,ui_old,vi_old = ufl.split(Ksol)
# ********* Electrical parameters ********* #
# ******** Save solution in Xdmf format ***
fileO = io.XDMFFile(msh.comm, "outE/ThermoElect.xdmf", "w")
fileO.write_mesh(msh)
# ********* Forcing terms and electrical stimulation ******* #
waveS0 = Function(Mh)
waveS0.interpolate(lambda x: stim_amp*(np.less(-6,x[0]))*(np.less(x[0],6))*(np.less(-6,x[1]))*(np.less(x[1],6))*(np.less(x[2],5)))
# Define the expression using FEniCSx
'''def WaveS0(amp, x):
return amp * (lt(-6, x[1]) and lt(-6, x[0]) and lt(x[2], 5))
# Create an instance of the expression
waveS0 = WaveS0(stim_amp,x)'''
# ***** Linear weak form ****
KL = ......
# initial condition
#Ksol.x.array[:] = 0.0
# assign previous solution
#Ksol0.x.array[:] = Ksol.x.array
bilinear_form = fem.form(KL)
LHS = assemble_matrix(bilinear_form)
LHS.assemble()
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(LHS)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
while (t <= Tfinal):
print("t=%.2f" % t)
KR = .......
#RHS = rhs(KR)
#LHS = lhs(KL)
#LHS1 = assemble_matrix(KL)
linear_form = fem.form(KR)
RHS = create_vector(linear_form)
# Update the right hand side reusing the initial vector
with RHS.localForm() as loc_b:
loc_b.set(0)
assemble_vector(RHS, linear_form)
# solver optimisation
# Solve linear problem
solver.solve(RHS, Ksol.vector)
ul,vl,ui,vi = Ksol.split()
Ksol.x.scatter_forward()
Ksol0.x.array[:] = Ksol.x.array
ul.name = "ul"; ui.name = "ui"; vl.name = "vl"; vi.name = "vi"
#fileO.write_function(Ksol, t)
if (inc % freqSave == 0):
fileO.write_function(ul, t)
fileO.write_function(vl, t)
fileO.write_function(vi, t)
fileO.write_function(ui, t)
t += dt;inc += 1
fileO.close()