How to output the data results on the boundary

Legacy DOLFIN
1

Hello everyone, I have solved the temperature and heat flux density distributions in a 2D domain using the following code. Suppose I want to output the heat flux density values of all nodes on a specific boundary (assuming the boundary’s index is 11) to a txt or excel file. How can I achieve this? I hope to obtain a general solution, and then I can use the same method to output the temperature on a specific boundary, and so on. Thank you!

# heat transfer

from dolfin import *
import matplotlib.pyplot as plt
import material_properties

mesh = Mesh('test.xml')      
boundary = MeshFunction("size_t", mesh, "test_facet_region.xml")
subdomains = MeshFunction("size_t", mesh, "test_physical_region.xml")

VT = FunctionSpace(mesh, "CG", 1)
T_ = TestFunction(VT)
dT = TrialFunction(VT)
Delta_T = Function(VT, name="Temperature increase")
ds = Measure('ds', domain=mesh, subdomain_data=boundary)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

q_gap0 = -0.035135     
g_m = Constant(q_gap0/material_properties.lambda_m)

aT_m = dot(grad(dT), grad(T_)) * dx(15)
LT_m = - g_m * T_ * ds(11)
bcT = [DirichletBC(VT, Constant(950.15), boundary, 10)]
solve(aT_m == LT_m, Delta_T,bcT)

submesh_m_T = SubMesh(mesh, subdomains, 15)
VT_m = FunctionSpace(submesh_m_T, 'P', 1)
Delta_T_m = interpolate(Delta_T, VT_m)      
vtkfile = File('m/T.pvd')
vtkfile << Delta_T_m

# heat flux field
Q = VectorFunctionSpace(mesh, 'P', 1)
w = TrialFunction(Q)
v = TestFunction(Q)
aQ = inner(w, v)*dx(15)
LQ = inner(-material_properties.lambda_m*grad(Delta_T), v)*dx(15)
flux = Function(Q)
solve(aQ == LQ, flux)

submesh_m_Q = SubMesh(mesh, subdomains, 15)
VQ_m = VectorFunctionSpace(submesh_m_Q, 'P', 1)
Q_m = interpolate(flux, VQ_m)
vtkfile = File('m/q.pvd')
vtkfile << Q_m
                       # Visualize flux field of the monolith
2

I would start with considering:

which you can restrict the facets in question with:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(2, 3)
V = FunctionSpace(mesh, "CG", 1)
u = project(Expression("x[0]+3*x[1]", degree=1), V)


class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)


mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
left = Left()
left.mark(mf, 1)
facets = mf.where_equal(1)


send_coordinates = []
send_values = []
send_global_idx = []

global_vertices = mesh.topology().global_indices(0)
vertex_values = u.compute_vertex_values()
coordinates = mesh.coordinates()

for facet_index in facets:
    facet = Facet(mesh, facet_index)
    vertices = facet.entities(0)
    for vertex in vertices:
        send_global_idx.append(global_vertices[vertex])
        send_coordinates.append(coordinates[vertex])
        send_values.append(vertex_values[vertex])
send_coordinates = np.asarray(
    send_coordinates).reshape(-1, coordinates.shape[1])

gathered_coordinates = mesh.mpi_comm().gather(send_coordinates, root=0)
gathered_values = mesh.mpi_comm().gather(send_values, root=0)
global_indices = mesh.mpi_comm().gather(send_global_idx, root=0)
if mesh.mpi_comm().rank == 0:
    unique_vertices = np.sort(np.unique(np.hstack(global_indices)))
    facet_values = np.zeros(len(unique_vertices))
    facet_coordinates = np.zeros((len(unique_vertices), coordinates.shape[1]))
    c = np.vstack(gathered_coordinates)
    v = np.hstack(gathered_values)
    i = np.hstack(global_indices)
    for coord, value, index in zip(c, v, i):
        local_idx = np.flatnonzero(unique_vertices == index)[0]
        facet_values[local_idx] = value
        facet_coordinates[local_idx] = coord
    for (coord, value) in zip(facet_coordinates, facet_values):
        print(coord, value)

yielding

 mpirun -n 4 python3 script.py
[1. 0.] 0.9999999999999999
[1.         0.33333333] 2.0000000000000004
[1.         0.66666667] 3.000000000000001
[1. 1.] 4.0000000000000036

You can easily store this to a text file using the normal file handler in Python

3

Thank you for your timely reply, dokken. I will study this method carefully

4

Dear dokken，I have another question about this. Suppose I obtain the heat flux density results on a 2D region boundary using your method. I want to use these results as the second type of boundary condition on a boundary of another 2D region, and then solve for the temperature field. How should I proceed? (The grids on these two boundaries are matched)

5

See for instance: How to impose a pressure distribution from a data file as a boundary condition (loading)? - #5 by dokken

6

Ok，I am learning the example that you have mentioned. I have a question: for two matching boundaries (in different domains) , I extract the computed results from one boundary and use them as the boundary conditions for the other boundary (assuming these results represent heat flux density); in fact, I simply assign the discretized data points at the mesh vertices to the corresponding vertices on the other boundary. When computing in the other domain, will FEniCS automatically interpolate these discrete data points into continuous boundary conditions before performing the computation? This is my understanding of the process, but I’m not sure if this is accurate. Can you help explain this to me?