Hello,
I’m having problems evaluating an integral with Dolfinx and I wonder if anyone could give me some tips. The integral is
in which u and T are the displacement and temperature of a body, obtained as the solution of a coupled thermoelastic problem. As a sanity check, I’m considering the case in which T = 0, so that the second and third term on the integral are expected to be zero. When I use assemble_scalar to evaluate each term separately, I get the expected result. However, when I assemble the whole expression, I get a different result. My question is why is this occurring.
I should also mention that I have two subdomains, which are identified from a level set function with UFL conditional, as discussed here. This is why the integral on D is obtained as the sum of the integral in each subdomain.
Here is the code which reproduces the error. It seems that the problem is in the last term, since when I replace T^2 with T, the error disappears.
Thanks in advance!
from dolfinx import fem
from ufl import FacetNormal, Measure, conditional, le, gt, dx, Identity, nabla_grad, div, inner, grad, TrialFunction, TestFunction, VectorElement, MixedElement, FiniteElement, split
from mpi4py import MPI
from dolfinx.mesh import (create_rectangle, CellType, locate_entities, meshtags)
import dolfinx.cpp as _cpp
from petsc4py.PETSc import ScalarType
import numpy as np
def tag_boundaries(msh, subdomains, dim):
indices, markers = [], []
for (marker, locator) in subdomains:
facets = locate_entities(msh, dim, locator)
indices.append(facets)
markers.append(np.full_like(facets, marker))
indices = np.hstack(indices).astype(np.int32)
markers = np.hstack(markers).astype(np.int32)
sorted_indices = np.argsort(indices)
return meshtags(msh, dim, indices[sorted_indices], markers[sorted_indices])
def apply_boundary_conditions(msh, bc_list):
fdim = msh.topology.dim - 1
bcs = []
for (facets, g, subspace, collapsed_subspace) in bc_list:
w_bc = fem.Function(collapsed_subspace)
w_bc.interpolate(lambda x: g(x))
bc_dofs = fem.locate_dofs_topological((subspace, collapsed_subspace), fdim, facets)
bcs.append(fem.dirichletbc(w_bc, bc_dofs, subspace))
return bcs
#---------------- Material properties ----------------------------------
E = 199.5E9; nu = 0.3;
mu, lbd= E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
alpha = 15.4e-6; kappa = 90; beta = alpha*E/(1-2*nu)
eps = 1e-3
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
def sigma(u):
return 2*mu*epsilon(u) + lbd*div(u)*Identity(cdim)
#------------------ Domain geometry -----------------------------------------------------
radius = 0.25; tol = 1e-8
msh = create_rectangle(comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (1.0, 1.0)), n=(200, 100),
cell_type=CellType.triangle, diagonal = _cpp.mesh.DiagonalType.crossed)
msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
comm = MPI.COMM_WORLD
cdim, fdim = msh.topology.dim, msh.topology.dim - 1
dx = Measure("dx", domain = msh); n = FacetNormal(msh)
#---------------- Material Applied loads ----------------------------------
Q = fem.Constant(msh, ScalarType( (0.0) ) )
F = fem.Constant(msh, ScalarType( (0.0, -2e7) ) )
#---------------- Function Spaces -------------------------------------------
degree_u = 1; degree_T = 1
el_u = VectorElement("CG", msh.ufl_cell(), degree_u, dim=2)
el_T = FiniteElement("CG", msh.ufl_cell(), degree_T)
el_m = MixedElement([el_u , el_T])
W = fem.FunctionSpace(msh, el_m)
U, _ = W.sub(0).collapse(); V, _ = W.sub(1).collapse()
# ------ Identifying the two-phases/subdomains from level-set ---------------
def phi(x):
return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius
phi_ = fem.Function(V)
phi_.interpolate(phi)
interior = conditional(le(phi_, tol), 1, 0); exterior = conditional(gt(phi_, tol), 1, 0)
# ------- Dirichlet boundary conditions -------------------------------------
u_dir = lambda x: (np.zeros_like(x[0]), np.zeros_like(x[0]))
T_dir = lambda x: (np.zeros_like(x[0]))
left_facets = locate_entities(msh, fdim, lambda x: np.isclose(x[0], 0.0))
right_facets = locate_entities(msh, fdim, lambda x: np.isclose(x[0], 1.0))
bc_list = [[np.hstack([left_facets, right_facets]), u_dir, W.sub(0), U],\
[np.hstack([left_facets, right_facets]), T_dir , W.sub(1), V]]
bcs = apply_boundary_conditions(msh, bc_list)
# ------- Neumann boundary conditions -------------------------------------
def NeumBd(x): return np.isclose(x[1], 0.0)
boundaries = [(1, lambda x: NeumBd(x))]
msh.topology.create_connectivity(msh.topology.dim - 1,msh.topology.dim)
facet_tags = tag_boundaries(msh, boundaries, fdim)
ds = Measure("ds", domain=msh, subdomain_data=facet_tags)
# ---- Solving thermoelastic problem -----------------------------------------
TrialF = TrialFunction(W); TestF = TestFunction(W)
(u, T) = split(TrialF); (v, r) = split(TestF)
a_state = inner(sigma(u),epsilon(v)) - beta * T * div(v) + kappa * inner(grad(T), grad(r))
lhs_state = a_state * interior * dx + eps * a_state * exterior * dx
rhs_state = inner(Q,r) * dx + inner(F,v) * ds(1)
petsc_opts={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
problem_state = fem.petsc.LinearProblem(lhs_state, rhs_state, bcs = bcs, petsc_options=petsc_opts)
wh = problem_state.solve()
(uh,Th) = wh.split(); uh = uh; Th = Th
print('\n--------- Assembling each term --------------------\n')
S1 = (1/2)*inner(sigma(uh),epsilon(uh))
term1 = fem.assemble_scalar(fem.form(S1*interior*dx + eps*S1*exterior*dx))
S1 = inner(beta, Th * div(uh))
term2 = fem.assemble_scalar(fem.form(S1*interior*dx + eps*S1*exterior*dx))
S1 = (1/2)* alpha * beta * Th**2 * 3.0
term3 = fem.assemble_scalar(fem.form(S1*interior*dx + eps*S1*exterior*dx))
print('term1 : '+str(term1))
print('term2 : '+str(term2))
print('term3 : '+str(term3)+'\n')
print('term1 - term2 : '+str(term1 - term2))
print('term1 - term2 + term3: '+str(term1 - term2 + term3))
print('\n--------- Assembling the whole expression ----------\n')
S1 = (1/2)*inner(sigma(uh),epsilon(uh))
compliance = fem.assemble_scalar(fem.form(S1*interior*dx + eps*S1*exterior*dx))
S1 = (1/2)*inner(sigma(uh),epsilon(uh)) - inner(beta, Th * div(uh))
mech_energy = fem.assemble_scalar(fem.form(S1*interior*dx + eps*S1*exterior*dx))
S1 = (1/2)*inner(sigma(uh),epsilon(uh)) - inner(beta, Th * div(uh)) + (1/2)* alpha * beta * Th**2 * 3.0
total_energy = fem.assemble_scalar(fem.form(S1*interior*dx + eps*S1*exterior*dx))
print('term1 : ' +str(compliance))
print('term1 - term2 : ' +str(mech_energy))
print('term1 - term2 + term3: ' +str(total_energy))``