Sure, as for the pde I’m trying to simulate, is the same I stated in this post.
My code is the following:
import numpy as np
import pyvista
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
div, dot, dx, grad, inner, lhs, rhs)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 3))
import ufl
import dolfinx
from petsc4py.PETSc import ScalarType
####
n_temp = 1.4
reff= -1.44*n_temp**-2+0.71*n_temp**-1+0.668+0.0636*n_temp
temp = (1+reff)/(1-reff)
A = Constant(mesh, ScalarType(temp))
####
x = ufl.SpatialCoordinate(mesh)
u_ex = lambda x: x[0]**3*x[1]+x[0]**2*x[1]+2*x[1]**2+10*x[1]
mu_s = Constant(mesh, ScalarType(1.22))
D = Constant(mesh, ScalarType(1/(3*1.22)))
mu_a = fem.Function(V)
mu_a.interpolate(lambda x: 0.036 + 0.082*(np.square(x[0]-0.5)+np.square(x[1]-0.3)<=np.square(0.2)))
n = ufl.FacetNormal(mesh)
g = -ufl.dot(D*ufl.grad(u_ex(x)),n)
S = fem.Constant(domain, ScalarType(0.0))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = D * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + mu_a * ufl.inner(u,v) * ufl.dx - ufl.inner(S, v) * ufl.dx
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = np.array(facet_tag.indices[facet_tag.values == marker])
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
boundary_conditions = [BoundaryCondition("Dirichlet", 1, lambda x: np.zeros(x[0].size)),
BoundaryCondition("Dirichlet", 2, lambda x: np.zeros(x[0].size)),
BoundaryCondition("Robin", 3, (1/(2*A), Constant(mesh, 0.0))),]
#BoundaryCondition("Neumann", 4, g)]
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
uh.name = 'uh'
#Data extraction
W = fem.VectorFunctionSpace(domain, ("DG", 0))
grad = ufl.grad(uh)
grad_expression = fem.Expression(grad, W.element.interpolation_points)
grad_h = fem.Function(W)
grad_h.interpolate(grad_expression)
outfile1 = open("dofs.csv", "w")
coords = V.tabulate_dof_coordinates().reshape((-1,gdim))
for x in coords:
print(x[0],",",x[1],file = outfile1)
##---------------------
outfile2 = open("solution.csv", "w")
sol = uh.x.array
for x in sol:
print(x, file = outfile2)
##------------------------
outfile3 = open("grad_solution.csv", "w")
???
##------------------------
outfile4 = open("mu_a.csv", "w")
data_mua = mu_a.x.array
for x in data_mua:
print(x, file = outfile4)
So what I’m trying to do is extract the dofs, the solution at the dofs, the gradient of the solution at the dofs and the constant mu_a at the dofs