jorgensd / jorgensd/dolfinx_mpc

"Non-square" MPC giving different sparsity pattern than square problem.

Open
#215 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
Python
Stars
62
Forks
22
Avg merge
1d 4h
Merged PRs (30d)
4

Description

```python
from mpi4py import MPI

import numpy as np
from dolfinx import fem, mesh
from basix.ufl import element
import dolfinx_mpc
import ufl

cell_type = mesh.CellType.hexahedron # BUG

N = 3
domain = mesh.create_unit_cube(MPI.COMM_WORLD, N, 2, N, cell_type=cell_type)

V = fem.functionspace(domain, element("Lagrange", domain.basix_cell(), 1))
Q = V.clone()

def periodic_boundary(x):
return np.logical_or(np.isclose(x[0], 1), np.isclose(x[2], 1))

def periodic_map(x):
out = x.copy()
out[0][np.isclose(x[0], 1).nonzero()] -= 1
out[2][np.isclose(x[2], 1).nonzero()] -= 1
return out

mpc_u = dolfinx_mpc.MultiPointConstraint(V)
mpc_u.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_map, [])
mpc_u.finalize()

mpc_p = dolfinx_mpc.MultiPointConstraint(Q)
mpc_p.create_periodic_constraint_geometrical(Q, periodic_boundary, periodic_map, [])
mpc_p.finalize()

# Stokes weak form
u, p = ufl.TrialFunction(V), ufl.TrialFunction(Q)
v, q = ufl.TestFunction(V), ufl.TestFunction(Q)

# Weak statementmpc_u,
a01 = v * p * ufl.dx
form_c = fem.form(a01)
mpcs = [mpc_u, mpc_p]

p0 = dolfinx_mpc.create_sparsity_pattern(form_c, mpcs)
p0.finalize()

p1 = dolfinx_mpc.create_sparsity_pattern(form_c, [mpc_u, mpc_u])
p1.finalize()
print(p0.num_nonzeros, p1.num_nonzeros)
breakpoint()
```
This in turn probably leads to an out of range error in assembly.

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.