Wonderful! Thanks for the help. I’d love to see this work!
EDIT: Did you manage to do the version change?
Ive not had time to do this yet.
Hello again,
Is this going to happen?
Hi @BillS
It’s not high on my list of priorities, as I’ve got more critical issues than making a backward compatible release of DOLFINx_MPC.
Best,
Jørgen
It’s OK. Sorry to be a bother. It was a little side project that I had hoped to get to work with as little headache as possible. I will find another solution (at least until maybe Dolfinx v0.6.0 appears in the Ubuntu PPA repository).
Cheers,
Bill
With the mid-August slowdown, I am revisiting this problem. The latest updates of my Ubuntu system included DolfinX version 0.6.0. I updated DolfinX-mpc from the PPA, but is did not work because it seems only the real version was available. So I uninstalled that and downloaded the MPC source version 0.6.0 and went through the compilation and installation process
cmake -G Ninja -DCMAKE_BUILD_TYPE=Release -B build-dir cpp/
ninja -j3 install -C build-dir
python3 -m pip install python/. --upgrade
which proceeded without error. I added the following environment variables to the BASH profile
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
export PATH=/usr/local/bin:$PATH
export PKG_CONFIG_PATH=/usr/local/lib/pkgconfig:$PKG_CONFIG_PATH
On running
# Generate unit cell of helical structure
import sys
import numpy as np
from mpi4py import MPI as nMPI
from petsc4py import PETSc
import slepc4py as slepc
import meshio
import ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx import io
import dolfinx_mpc
import gmsh
l = 6.0 # Cell length
r = 6.0 # Cell radius
rh = 2.0 # Helix radius
rw = 0.2 # wire radius
npts = 100 # number of points of spline curve
eps = 1.0e-5
Gamma = 0 + 6.283j
BCweight = 100.0
comm = nMPI.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # process where GMSG runs
if mpiRank == modelRank:
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)
gmsh.model.add("Helix Cell")
c = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, l, r, 1, 2*np.pi)
dt = 2 * np.pi / npts
p = []
for n in range(npts+3):
theta = 2 * np.pi * (n-1) / npts
kk = gmsh.model.occ.addCircle(rh * np.cos(theta), rh * np.sin(theta), (n -1) * l / npts, rw)
p.append(gmsh.model.occ.addCurveLoop([kk], tag=-1))
print(p)
s = gmsh.model.occ.addThruSections(p, tag=-1)
print(s)
v = gmsh.model.occ.cut([(3,c)], s, tag=100, removeObject=True, removeTool=True)
print("v={}".format(v))
for n in range(npts+3):
gmsh.model.occ.remove([(1,p[n])], recursive=True)
gmsh.model.occ.synchronize()
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
print(pt)
# Boundary conditions
Bot = []
Top = []
PEC = []
for bb in pt:
CoM = gmsh.model.occ.getCenterOfMass(bb[0], bb[1])
print(CoM)
if(np.abs(CoM[2]) < eps):
Bot.append(bb[1])
elif(np.abs(CoM[2]-l) < eps):
Top.append(bb[1])
else:
PEC.append(bb[1])
print(Bot, Top, PEC)
gmsh.model.mesh.setPeriodic(2, Top, Bot, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, l, 0, 0, 0, 1])
gmsh.model.addPhysicalGroup(3, [100], 1) # Helix volume
gmsh.model.setPhysicalName(3, 1, "HelixVolume")
gmsh.model.addPhysicalGroup(2, Top, 1) # Top BC (slave periodic BC)
gmsh.model.addPhysicalGroup(2, Bot, 2) # Bottom BC (master periodic BC)
gmsh.model.addPhysicalGroup(2, PEC, 3) # PEC boundary
gmsh.model.setPhysicalName(2, 1, "TopPBC")
gmsh.model.setPhysicalName(2, 2, "BotPBC")
gmsh.model.setPhysicalName(2, 3, "PEC")
pt = gmsh.model.getEntitiesInBoundingBox(-eps-rh-rw, -rh-rw-eps, -eps, rh+rw+eps, rh+rw+eps, eps)
gmsh.model.mesh.setSize(pt, 0.05)
print(pt)
pt = gmsh.model.getEntitiesInBoundingBox(-eps-rh-rw, -rh-rw-eps, l-eps, rh+rw+eps, rh+rw+eps, l+eps)
gmsh.model.mesh.setSize(pt, 0.05)
print(pt)
#gmsh.model.mesh.setSize([(0,108),(0,109)], 0.05)
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.001)
gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)
gmsh.option.setNumber('Mesh.Format', 1)
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 50)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.model.mesh.generate(3)
#gmsh.fltk.run()
# mesh = 3d topology
# ct = volume tags
# fm = surface tags
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if modelRank == mpiRank:
gmsh.finalize()
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
V = fem.VectorFunctionSpace(mesh, ("CG", 2))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Print out BC tags
with io.XDMFFile(mesh.comm, "BCs.xdmf", "w") as xx:
xx.write_mesh(mesh)
xx.write_meshtags(fm)
n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
# Periodic BC mapping
def PeriodicMapping(x):
out = np.zeros(x.shape)
out[0] = x[0]
out[1] = x[1]
out[2] = x[2] + l
return out
# Dirichlet BCs using weak invokation
Lb = BCweight * ufl.inner(ufl.cross(n, u), ufl.cross(n, v)) * ds(3)
bc=[]
# Create multipoint boundary. Master boundary has tag 2,
# slave has tag 1 - set on periodic constraint.
# (see GMSH physical group definitions above)
fbot = fm.find(2)
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V, fm, 1, PeriodicMapping, bc, 1.0)
mpc.finalize()
# Set up forms for eigenvalue problem Ax -k0^2*Bx = 0
az = ufl.as_vector((0, 0, 1)) # Z-directed unit vector
a = fem.form(ufl.inner(Gamma * ufl.cross(az, u) + ufl.curl(u), Gamma * ufl.cross(az, v) + ufl.curl(v)) * ufl.dx + Lb)
b = fem.form(ufl.inner(u, v) * ufl.dx)
num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
if mpiRank == modelRank:
print(f"Number of dofs global: {num_dofs_global}")
A = dolfinx_mpc.assemble_matrix(a, mpc, bcs=bc, diagval=1.0)
B = dolfinx_mpc.assemble_matrix(b, mpc, bcs=bc, diagval=1.0)
eigs = slepc.EPS().create(mesh.comm)
eigs.setOperators(A, B)
eigs.setProblemType(slepc.EPS.ProblemType.GNHEP) # General non-hermitian
eigs.setTolerances(1.0e-8)
eigs.setType(slepc.EPS.Type.KRYLOVSCHUR)
st = eigs.getST() # Invoke spetcral transformation
st.setType(slepc.ST.Type.SINVERT)
eigs.setWhichEigenpairs(slepc.EPS.Which.TARGET_REAL) # target real part of eigenvalue
eigs.setTarget(-1.0)
eigs.setDimensions(nev=5) # Number of eigenvalues
eigs.solve()
eigs.view()
eigs.errorView()
sys.exit(0)
I get an error that says dolfinx-mpc does not exist. What am I missing here?
Traceback (most recent call last):
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 13, in <module>
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 13, in <module>
import dolfinx_mpc
import dolfinx_mpc
ModuleNotFoundError: No module named 'dolfinx_mpc'
ModuleNotFoundError: No module named 'dolfinx_mpc'
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 13, in <module>
import dolfinx_mpc
ModuleNotFoundError: No module named 'dolfinx_mpc'
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 13, in <module>
import dolfinx_mpc
ModuleNotFoundError: No module named 'dolfinx_mpc'
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 13, in <module>
import dolfinx_mpc
ModuleNotFoundError: No module named 'dolfinx_mpc'
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 13, in <module>
import dolfinx_mpc
ModuleNotFoundError: No module named 'dolfinx_mpc'
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 13, in <module>
import dolfinx_mpc
ModuleNotFoundError: No module named 'dolfinx_mpc'
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 13, in <module>
import dolfinx_mpc
ModuleNotFoundError: No module named 'dolfinx_mpc'
--------------------------------------------------------------------------
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:
Process name: [[16316,1],5]
Exit code: 1
--------------------------------------------------------------------------
Start by checking that
python3 -m pip list
shows dolfinx_mpc.
Next you would need to post the output of the installation commands.
Pip did not show dolfinx_mpc. In fact, I found out that I had downloaded the wrong file from Github (the Dolfinx blob instead of dolfinx-mpc…boy do I feel silly). In fact, at the installation stage, it did not actually install anything, since dolfinx 0.6.0 is already installed.
Pip list now shows dolfinx-mpc and the program runs up until mpc.finalize(), where I get the following cryptic error:
Traceback (most recent call last):
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 138, in <module>
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 138, in <module>
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 138, in <module>
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 138, in <module>
mpc.finalize()
mpc.finalize()
File "/usr/lib/python3/dist-packages/dolfinx_mpc/multipointconstraint.py", line 97, in finalize
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 138, in <module>
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 138, in <module>
mpc.finalize()
File "/usr/lib/python3/dist-packages/dolfinx_mpc/multipointconstraint.py", line 97, in finalize
mpc.finalize()
File "/usr/lib/python3/dist-packages/dolfinx_mpc/multipointconstraint.py", line 97, in finalize
self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
File "/home/bill/.local/lib/python3.10/site-packages/dolfinx/fem/function.py", line 469, in __init__
self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
File "/home/bill/.local/lib/python3.10/site-packages/dolfinx/fem/function.py", line 469, in __init__
mpc.finalize()
File "/usr/lib/python3/dist-packages/dolfinx_mpc/multipointconstraint.py", line 97, in finalize
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 138, in <module>
self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
File "/home/bill/.local/lib/python3.10/site-packages/dolfinx/fem/function.py", line 469, in __init__
self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
File "/home/bill/.local/lib/python3.10/site-packages/dolfinx/fem/function.py", line 469, in __init__
mpc.finalize()
File "/usr/lib/python3/dist-packages/dolfinx_mpc/multipointconstraint.py", line 97, in finalize
File "/usr/lib/python3/dist-packages/dolfinx_mpc/multipointconstraint.py", line 97, in finalize
self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
File "/home/bill/.local/lib/python3.10/site-packages/dolfinx/fem/function.py", line 469, in __init__
mpc.finalize()
File "/usr/lib/python3/dist-packages/dolfinx_mpc/multipointconstraint.py", line 97, in finalize
self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
File "/home/bill/.local/lib/python3.10/site-packages/dolfinx/fem/function.py", line 469, in __init__
self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
File "/home/bill/.local/lib/python3.10/site-packages/dolfinx/fem/function.py", line 469, in __init__
assert mesh is None
AssertionError
assert mesh is None
AssertionError
assert mesh is None
AssertionError
assert mesh is None
AssertionError
assert mesh is None
AssertionError
assert mesh is None
AssertionError
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 138, in <module>
mpc.finalize()
File "/usr/lib/python3/dist-packages/dolfinx_mpc/multipointconstraint.py", line 97, in finalize
self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
File "/home/bill/.local/lib/python3.10/site-packages/dolfinx/fem/function.py", line 469, in __init__
assert mesh is None
AssertionError
assert mesh is None
AssertionError
--------------------------------------------------------------------------
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
Could you show the output of python3 -m pip list, i.e. show what versions of DOLFINx and DOLFINx_mpc you are using?
Here you go!
Cheers
Package Version
------------------------ ------------------------------------
alabaster 0.7.12
appdirs 1.4.4
apturl 0.5.2
attrs 21.2.0
autobahn 21.11.1
Automat 20.2.0
Babel 2.8.0
base58 1.0.3
bcrypt 3.2.0
beautifulsoup4 4.10.0
beniget 0.4.1
blinker 1.4
Brlapi 0.8.3
Brotli 1.0.9
cbor 1.0.0
certifi 2020.6.20
cffi 1.15.0
cftime 1.5.2
chardet 4.0.0
click 8.0.3
colorama 0.4.4
command-not-found 0.3
commonmark 0.9.1
constantly 15.1.0
cryptography 3.4.8
cupshelpers 1.0
cycler 0.11.0
Cython 0.29.28
dbus-python 1.2.18
decorator 4.4.2
defer 1.0.6
distro 1.7.0
distro-info 1.1+ubuntu0.1
docutils 0.17.1
dolfinx-mpc 0.6.1
duplicity 0.8.21
ecdsa 0.18.0b1
fasteners 0.14.1
fenics-basix 0.6.0
fenics-dijitso 2019.2.0.dev0
fenics-dolfin 2019.2.0.dev0
fenics-dolfinx 0.6.0
fenics-ffc 2019.2.0.dev0
fenics-ffcx 0.6.0
fenics-fiat 2019.2.0.dev0
fenics-ufl 2023.1.1.post0
fenics-ufl-legacy 2022.3.0
flatbuffers 1.12.1-git20200711.33e2d80-dfsg1-0.6
fonttools 4.29.1
fs 2.4.12
future 0.18.2
gast 0.5.2
gpg 1.16.0
h5py 3.6.0
h5py.-debian-h5py-serial 3.6.0
html5lib 1.1
httplib2 0.20.2
hyperlink 21.0.0
idna 3.3
imageio 2.23.0
imagesize 1.3.0
importlib-metadata 4.6.4
incremental 21.3.0
jeepney 0.7.1
Jinja2 3.0.3
keyring 23.5.0
kiwisolver 1.3.2
language-selector 0.1
launchpadlib 1.10.16
lazr.restfulclient 0.14.4
lazr.uri 1.0.6
llvmlite 0.38.0
lockfile 0.12.2
louis 3.20.0
lxml 4.8.0
lz4 3.1.3+dfsg
macaroonbakery 1.3.1
Mako 1.1.3
MarkupSafe 2.0.1
matplotlib 3.5.1
meshio 5.3.0
mnemonic 0.19
monotonic 1.6
more-itertools 8.10.0
mpi4py 3.1.3
mpmath 0.0.0
msgpack 1.0.3
mshr 2019.2.0.dev0
netCDF4 1.5.8
netifaces 0.11.0
numba 0.55.1
numpy 1.21.5
oauthlib 3.2.0
olefile 0.46
packaging 21.3
paramiko 2.9.3
passlib 1.7.4
petsc4py 3.15.1
pexpect 4.8.0
Pillow 9.0.1
pip 22.0.2
pkgconfig 1.5.5
ply 3.11
pooch 1.6.0
protobuf 3.12.4
ptyprocess 0.7.0
py-ubjson 0.16.1
pyasn1 0.4.8
pyasn1-modules 0.2.1
pybind11 2.9.1
pycairo 1.20.1
pycparser 2.21
pycups 2.0.1
Pygments 2.11.2
PyGObject 3.42.1
PyHamcrest 2.0.2
PyJWT 2.3.0
pymacaroons 0.13.0
PyNaCl 1.5.0
pyOpenSSL 21.0.0
pyparsing 2.4.7
pypng 0.0.20
PyQRCode 1.2.1
pyRFC3339 1.1
python-apt 2.4.0+ubuntu2
python-dateutil 2.8.1
python-debian 0.1.43+ubuntu1.1
python-snappy 0.5.3
pythran 0.10.0
PyTrie 0.4.0
pytz 2022.1
pyvista 0.37.0
pyxdg 0.27
PyYAML 5.4.1
reportlab 3.6.8
requests 2.25.1
rich 11.2.0
roman 3.3
scipy 1.8.0
scooby 0.7.0
SecretStorage 3.3.1
service-identity 18.1.0
setuptools 59.6.0
six 1.16.0
slepc4py 3.15.1
snowballstemmer 2.2.0
sortedcontainers 2.1.0
soundconverter 4.0.3
soupsieve 2.3.1
Sphinx 4.3.2
ssh-import-id 5.11
sympy 1.9
systemd-python 234
Twisted 22.1.0
txaio 21.2.1
u-msgpack-python 2.3.0
ubuntu-advantage-tools 8001
ubuntu-drivers-common 0.0.0
ufoLib2 0.13.1
ufw 0.36.1
ujson 5.1.0
unattended-upgrades 0.1
unicodedata2 14.0.0
urllib3 1.26.5
usb-creator 0.3.7
vtk 9.2.4
wadllib 1.3.6
webencodings 0.5.1
wheel 0.37.1
wsaccel 0.6.3
xdg 5
xkit 0.0.0
zipp 1.0.0
zope.interface 5.4.0
I’m not able to reproduce your error with a slightly modified script:
docker run -ti -v $(pwd):/root/shared -w /root/shared --rm ghcr.io/jorgensd/dolfinx_mpc:v0.6.1.post1
and
import dolfinx_mpc
import gmsh
import numpy as np
import ufl
from dolfinx import fem, io
from dolfinx.io import gmshio
from mpi4py import MPI as nMPI
from slepc4py import SLEPc as slepc
l = 6.0 # Cell length
r = 6.0 # Cell radius
rh = 2.0 # Helix radius
rw = 0.2 # wire radius
npts = 100 # number of points of spline curve
eps = 1.0e-5
Gamma = 0 + 6.283j
BCweight = 100.0
comm = nMPI.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # process where GMSG runs
if mpiRank == modelRank:
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)
gmsh.model.add("Helix Cell")
c = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, l, r, 1, 2*np.pi)
dt = 2 * np.pi / npts
p = []
for n in range(npts+3):
theta = 2 * np.pi * (n-1) / npts
kk = gmsh.model.occ.addCircle(rh * np.cos(theta), rh * np.sin(theta), (n -1) * l / npts, rw)
p.append(gmsh.model.occ.addCurveLoop([kk], tag=-1))
print(p)
s = gmsh.model.occ.addThruSections(p, tag=-1)
print(s)
v = gmsh.model.occ.cut([(3,c)], s, tag=100, removeObject=True, removeTool=True)
print("v={}".format(v))
for n in range(npts+3):
gmsh.model.occ.remove([(1,p[n])], recursive=True)
gmsh.model.occ.synchronize()
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
print(pt)
# Boundary conditions
Bot = []
Top = []
PEC = []
for bb in pt:
CoM = gmsh.model.occ.getCenterOfMass(bb[0], bb[1])
print(CoM)
if(np.abs(CoM[2]) < eps):
Bot.append(bb[1])
elif(np.abs(CoM[2]-l) < eps):
Top.append(bb[1])
else:
PEC.append(bb[1])
print(Bot, Top, PEC)
gmsh.model.mesh.setPeriodic(2, Top, Bot, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, l, 0, 0, 0, 1])
gmsh.model.addPhysicalGroup(3, [100], 1) # Helix volume
gmsh.model.setPhysicalName(3, 1, "HelixVolume")
gmsh.model.addPhysicalGroup(2, Top, 1) # Top BC (slave periodic BC)
gmsh.model.addPhysicalGroup(2, Bot, 2) # Bottom BC (master periodic BC)
gmsh.model.addPhysicalGroup(2, PEC, 3) # PEC boundary
gmsh.model.setPhysicalName(2, 1, "TopPBC")
gmsh.model.setPhysicalName(2, 2, "BotPBC")
gmsh.model.setPhysicalName(2, 3, "PEC")
pt = gmsh.model.getEntitiesInBoundingBox(-eps-rh-rw, -rh-rw-eps, -eps, rh+rw+eps, rh+rw+eps, eps)
gmsh.model.mesh.setSize(pt, 0.05)
pt = gmsh.model.getEntitiesInBoundingBox(-eps-rh-rw, -rh-rw-eps, l-eps, rh+rw+eps, rh+rw+eps, l+eps)
gmsh.model.mesh.setSize(pt, 0.05)
#gmsh.model.mesh.setSize([(0,108),(0,109)], 0.05)
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.001)
gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)
gmsh.option.setNumber('Mesh.Format', 1)
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 50)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.model.mesh.generate(3)
comm.Barrier()
# mesh = 3d topology
# ct = volume tags
# fm = surface tags
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if modelRank == mpiRank:
gmsh.finalize()
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
V = fem.VectorFunctionSpace(mesh, ("CG", 2))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Print out BC tags
with io.XDMFFile(mesh.comm, "BCs.xdmf", "w") as xx:
xx.write_mesh(mesh)
xx.write_meshtags(fm)
n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
# Periodic BC mapping
def PeriodicMapping(x):
out = np.zeros(x.shape)
out[0] = x[0]
out[1] = x[1]
out[2] = x[2] + l
return out
# Dirichlet BCs using weak invokation
Lb = BCweight * ufl.inner(ufl.cross(n, u), ufl.cross(n, v)) * ds(3)
bc=[]
# Create multipoint boundary. Master boundary has tag 2,
# slave has tag 1 - set on periodic constraint.
# (see GMSH physical group definitions above)
fbot = fm.find(2)
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V, fm, 1, PeriodicMapping, bc, 1.0)
mpc.finalize()
# Set up forms for eigenvalue problem Ax -k0^2*Bx = 0
az = ufl.as_vector((0, 0, 1)) # Z-directed unit vector
a = fem.form(ufl.inner(Gamma * ufl.cross(az, u) + ufl.curl(u), Gamma * ufl.cross(az, v) + ufl.curl(v)) * ufl.dx + Lb)
b = fem.form(ufl.inner(u, v) * ufl.dx)
num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
if mpiRank == modelRank:
print(f"Number of dofs global: {num_dofs_global}")
A = dolfinx_mpc.assemble_matrix(a, mpc, bcs=bc, diagval=1.0)
B = dolfinx_mpc.assemble_matrix(b, mpc, bcs=bc, diagval=1.0)
eigs = slepc.EPS().create(mesh.comm)
eigs.setOperators(A, B)
eigs.setProblemType(slepc.EPS.ProblemType.GNHEP) # General non-hermitian
eigs.setTolerances(1.0e-8)
eigs.setType(slepc.EPS.Type.KRYLOVSCHUR)
st = eigs.getST() # Invoke spetcral transformation
st.setType(slepc.ST.Type.SINVERT)
eigs.setWhichEigenpairs(slepc.EPS.Which.TARGET_REAL) # target real part of eigenvalue
eigs.setTarget(-1.0)
eigs.setDimensions(nev=5) # Number of eigenvalues
eigs.solve()
eigs.view()
eigs.errorView()
or with the following conda environment:
name: mpc_test
channels:
- conda-forge
dependencies:
- fenics-dolfinx
- dolfinx_mpc
- petsc=*=complex*
- python-gmsh
It is certainly something to do with my installation. I will need to dig into it a bit - maybe do a complete purge and reinstall. Thanks for your help and have a great weekend!