Dear all, I am trying to solve a Poisson eqaution with DG spaces, in a region given by a square with a circle in it. The circle is meshed inside.
If I define my variational problem (here I am showing only one part of it) as
F_A = (fsp.u.dx(i) * fsp.nu_u.dx(i)) * dx + \
(fsp.f * fsp.nu_u) * dx + \
- facet_normal[i] * (fsp.u.dx(i)) * fsp.nu_u * ds
then it works, I get the exact solution.
But if I use
F_B = (fsp.u.dx(i) * fsp.nu_u.dx(i)) * dx + \
(fsp.f_shape * fsp.nu_u) * dx_shape + \
(fsp.f_square * fsp.nu_u) * dx_square + \
- facet_normal[i] * (fsp.u.dx(i)) * fsp.nu_u * ds
I get the wrong one. Note that dx, dx_shape, dx_square have been tested carefully by integrating a function on them, and they are correct.
I set the f_shape and f_square profiles with
msh.interpolate_dg(fsp.f, f_shape_expression(), sf, shape_id)
msh.interpolate_dg(fsp.f, f_square_expression(), sf, square_id)
msh.interpolate_dg(fsp.f_shape, f_shape_expression())
msh.interpolate_dg(fsp.f_square, f_square_expression())
where interpolate_dg is
def interpolate_dg(f, g, sf=None, region_id=None):
Q = f.function_space()
element = Q.ufl_element()
if (element.family() != 'Discontinuous Lagrange'):
# the method has been called on a field defined on a continuous function space -> throw an error and exit
print(f'{col.Fore.RED}Error: interpolate_dg has been called on a field with a continuous function space!! Stopping now.{col.Fore.RESET}')
sys.exit(1)
if (element.value_shape() != g.value_shape()):
# the value shape of Q and that of g differ -> check whether this is due to a 'convention' issue where scalars have been given a value shape of (1,) vs. ()
if ((((element.value_shape() == ()) and (g.value_shape() == (1,))) or ((element.value_shape() == (1,)) and (g.value_shape() == ()))) == False):
# the discrepancy was not due to a convention issue -> throw an error and exit
print(f'{col.Fore.RED}Error: value shapes are different!!\n\telement value shape = {element.value_shape()}\n\tg value shape= {g.value_shape()}\nStopping now.{col.Fore.RESET}')
sys.exit(1)
value_size = int(np.prod(element.value_shape())) if element.value_shape() else 1
mesh = Q.mesh()
'''
dof_coordinates stores the coordinates of the points where DOFs sit. Because the field 'f' defined on each DOF has value_size components, dof_coordinates is composed of blocks, where each block has 'value_size' entries, and blocks are all identical
For example, dof_coordinates is of the form ->
row 0: [x0, y0] ← this corresponds to f[0] at DOF point 0
row 1: [x0, y0] ← this corresponds to f[1] at DOF point 0
...
row value_size [x1, y1] ← this corresponds to f[0] at DOF point 1
row value_size+1 [x1, y1] ← this corresponds to f[1] at DOF point 1
...
'''
dof_coordinates = Q.tabulate_dof_coordinates()
'''
f_values contains the value of 'f' on the DOFs, and it has the same structure as 'dof_coordinates'
f_values =
entry 0: f[0] at DOF point 0
entry 1: f[1] at DOF point 0
...
entry value_size f[0] at DOF point 1
entry value_size + 1 f[1] at DOF point 1
...
'''
f_values = f.vector().get_local() # get a copy of field values
for cell in cells(mesh):
# run on all mesh cells
if region_id != None:
#region_id has been given when calling this method -> evaluate sf on the cell to obtain the tag of cell `cell`
# compute 'sf' on the cell; under consideration
cell_tag = sf[cell]
'''
cell_dofs contains the IDs of the DOFs that are contained into 'cell', it has the structure
[
id_f_0_on_DOF_0,
id_f_0_on_DOF_1,
...,
id_f_0_on_DOF_{n_nodes-1},
id_f_1_on_DOF_0,
id_f_1_on_DOF_1,
...,
id_f_1_on_DOF_{n_nodes-1},
...
]
where the pattern is repeated value_size times, i.e., one for each component of 'f', and n_nodes = [number of DOFs in the cell] / [value_size]. In other words
cell_dofs[j * n_nodes + i] = [index in f.values().get_local() corresponding to the j-th component of the field 'f' sitting on ith DOF in the cell 'cell']
'''
cell_dofs = Q.dofmap().cell_dofs(cell.index())
n_nodes = len(cell_dofs) // value_size
'''
remove the redundancy in cell_dofs and store the result in
cell_dofs_unique =
[
id_DOF_0,
id_DOF_1,
...,
id_DOF_{n_nodes-1}
]
'''
cell_dofs_unique = cell_dofs[:n_nodes]
if (region_id == None) or (cell_tag == region_id):
# if 'cell_tag' == 'id', then the cell under consideration belongs to the surface tagged with 'id' -> set the DOFs of 'f' relative to this cell according to 'g'
for i in range(len(cell_dofs_unique)):
# run over physical DOFs contained to 'cell' and print out the value of 'f' by specifying that those DOFs belong to region tagged with 'cell_tag' in a separate column of the csv output file. Note that, because the space of 'f' is discontinuous, here DOFs in 'cell' may belong to different mesh regions, and thus have different tags
# pad 'x' to three dimensions
dof_coordinate = dof_coordinates[cell_dofs_unique[i]]
for j in range(value_size):
f_values[cell_dofs[j * n_nodes + i]] = np.atleast_1d(g(dof_coordinate[:2]))[j]
f.vector().set_local(f_values)
f.vector().apply("insert")
Why simply replacing
(fsp.f * fsp.nu_u) * dx
with
(fsp.f_shape * fsp.nu_u) * dx_shape + (fsp.f_square * fsp.nu_u) * dx_square
changes the result?
I can provide more details about the code, if needed.