Open florian98765 opened 3 years ago
wence-
commented
3 years ago 1246 suggests that petsc is able to ready periodic meshes, but one needs to create a proper discontinuous coordinate field.
I think this is the right thing to fix if one wants a periodic topology.
In some circumstances (e.g. anti-periodic domains) one doesn't want a periodic topology. The way to do that is to put the periodicity into the globaltolocal maps. This could also be done.
pefarrell
commented
3 years ago I have implemented an approach like DOLFIN's in Firedrake. I identify matching master-slave dofs from user-specified data in a preprocessing step, then adjust the assembled equations using the post-assembly hooks. The approach is parallelised but could be cleaned up/generalised. (I couldn't use the built-in periodic constructors because I needed to enforce an anti-periodic condition, u|left = -1 * u|right).
Term is starting and I don't really have time to take this code forward. Perhaps I could share it with you, if you'd be willing to abstract the approach and make a PR?
medinaeder
commented
3 years ago This will be incredibly useful to me as well and I would be willing to give it a shot.
@pefarrell is this code similar to @florianwechsung firedrake-periodicity? Currently it is missing parallel, mg, and antiperiodic support.
florian98765
commented
3 years ago Sounds great. Maybe we can work together on a cleaned up and generalized version. How can you share the code?
Thanks a lot and best wishes Florian
florianwechsung
commented
3 years ago I think what @pefarrell is doing in his code is different from firedrake-periodicity. I create a periodic mesh object and then a dg coordinate field. From what Patrick describes, he builds the matrices on the original mesh, and then modifies them afterwards.
I don't think firedrake-periodicity can straightforwardly be made to support antiperiodic bcs, but parallel and mg should not be a huge amount of work. Essentially you just have to make sure that you perform the mesh modification before the mesh is distributed.
pefarrell
commented
3 years ago Great! Here it is.
The important parts are in def function_space(...) and def solver(...). In function_space, I identify the master-slave dof pairs in global numbering, using the data provided:
masterslave is a Function that is negative for slaves, positive for masterssign is a Function that says whether to be periodic or antiperiodic (negative for periodic)key is a Python function that takes in an index in local mixed numbering and returns the coordinates that vary over the face I'm trying to match up, as well as the vector-component for vector-valued spaces. In this example I have a 2D box and I want to match x = -1/2 with x = +1/2, so I return the y-component of the coordinate field.The setup code then creates various data structures for later use: a scatter from master to slave across processes, lists of matching slave and master dofs in global numbering, and a temporary PETSc Vec to use with the scatter.
Once all that is done, the modification of the assembled matrices is fairly easy. It happens in post_function_callback and post_jacobian_callback. We modify the equation for slave dofs to be X[slave] + sign[slave] * X[master] = 0, and modify the rows of the Jacobian accordingly.
This same process would need to happen on every level of a multigrid hierarchy. It also means that you can't use solvers that rely purely on symbolic information (e.g. PCPATCH, although the new ASM mechanism would replace it).
Lastly we should probably change the nomenclature away from master/slave, cf https://gizmodo.com/github-to-remove-master-and-slave-coding-terms-widely-s-1844041329 .
The problem I'm solving here is difficult and it takes 218 SNES iterations on my machine to converge. Sorry about this! I didn't have time to set up a simpler problem.
# -*- coding: utf-8 -*-
from firedrake import *
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
class SmecticProblem(object):
def mesh(self, comm):
baseN = 30
base = UnitSquareMesh(3*baseN, baseN, comm=comm, quadrilateral=True)
base.coordinates.dat.data[:,0] -= 0.5
return base
def function_space(self, mesh):
U = FunctionSpace(mesh, "CG", 3)
V = VectorFunctionSpace(mesh, "CG", 2, dim=2)
Z = MixedFunctionSpace([U, V])
print("Z.dim(): %s %s" % (Z.dim(), [Z.sub(i).dim() for i in range(2)]))
# Data for enforcing periodic boundary conditions
comm = mesh.comm
# Function to track relation: negative to indicate slave, positive to indicate master
masterslave = Function(Z)
DirichletBC(Z.sub(0), -1, 1).apply(masterslave)
DirichletBC(Z.sub(1), Constant((-1, -1)), 1).apply(masterslave)
DirichletBC(Z.sub(0), +1, 2).apply(masterslave)
DirichletBC(Z.sub(1), Constant((+1, +1)), 2).apply(masterslave)
# Sign of the linear combination to take: if the value at slave dof is
# alpha, equation to be enforced there is X[slave] + alpha * X[master] = 0
sign = Function(Z)
DirichletBC(Z.sub(0), -1, 1).apply(sign) # periodic in density
DirichletBC(Z.sub(1).sub(0), -1, 1).apply(sign) # periodic in q1
DirichletBC(Z.sub(1).sub(1), +1, 1).apply(sign) # antiperiodic in q2
# Key function. This takes in a _local mixed index_
# and returns (relevant coordinates) + (index in bs)
# This is used to match up master and slave dofs.
num_density = Z._ises[0].getLocalSize()
bs = V.value_size
coordU = VectorFunctionSpace(mesh, "CG", 3, dim=2)
coordV = VectorFunctionSpace(mesh, "CG", 2, dim=2)
X = SpatialCoordinate(mesh)
XU = interpolate(X, coordU)
XV = interpolate(X, coordV)
def key(i):
if i < num_density:
# i refers to a density dof
return (XU.dat.data_ro[i, 1], -1)
else:
# i refers to a director dof
return (XV.dat.data_ro[(i - num_density) // bs, 1], (i - num_density) % bs)
# Do allgather to swap data
sendobj = {"master": [], "slave": []}
(lo, hi) = Z.dof_dset.layout_vec.getOwnershipRange()
with masterslave.dat.vec_ro as msv:
msarray = msv.array
for local_mixed_idx in range(hi - lo):
ms = msarray[local_mixed_idx]
if ms == 0:
continue
global_mixed_idx = lo + local_mixed_idx
keydata = key(local_mixed_idx)
if ms < 0:
sendobj["slave"].append((global_mixed_idx, keydata))
else:
sendobj["master"].append((global_mixed_idx, keydata))
recvobj = comm.allgather(sendobj)
unifyobj = {}
for entry in recvobj:
for label in entry:
unifyobj[label] = unifyobj.get(label, []) + entry[label]
for label in unifyobj:
unifyobj[label] = sorted(unifyobj[label], key=lambda x: x[1])
for (x, y) in zip(unifyobj["slave"], unifyobj["master"]):
#print("%s: global dofs %s %s: coordinates %s %s" % (comm.rank, x[0], y[0], x[1], y[1]), flush=True)
assert x[1] == y[1]
all_slave_global_mixed_indices = [x[0] for x in unifyobj["slave"]]
all_master_global_mixed_indices = [x[0] for x in unifyobj["master"]]
matching_indices = [(slave, master) for (slave, master) in zip(all_slave_global_mixed_indices, all_master_global_mixed_indices) if lo <= slave < hi]
if len(matching_indices) > 0:
my_slave_global_mixed_indices, my_master_global_mixed_indices = zip(*matching_indices)
else:
my_slave_global_mixed_indices, my_master_global_mixed_indices = [], []
# Now set up the scatter from master -> slave
Xscat = Z.dof_dset.layout_vec.duplicate()
slaveis = PETSc.IS().createGeneral(my_slave_global_mixed_indices, comm=comm)
masteris = PETSc.IS().createGeneral(my_master_global_mixed_indices, comm=comm)
Xtmp = Xscat.duplicate()
scatter = PETSc.Scatter().create(Xtmp, masteris, Xscat, slaveis)
# Store stuff necessary for assembly
self.sign = sign
self.my_slave = my_slave_global_mixed_indices
self.my_master = my_master_global_mixed_indices
self.scatter = scatter
self.Xscat = Xscat
return Z
def energy(self, z, params):
q = params[0]
W = params[1]
r = params[2]
a = Constant(-5*2)
b = Constant(0)
c = Constant(5*2)
B = Constant(1e-5)
K = Constant(0.3)
l = Constant(30)
s = FacetNormal(z.function_space().mesh())
(u, d) = split(z)
Q = as_tensor([[d[0], d[1]],
[d[1], -d[0]]])
Q_bottom = as_tensor([[1/2, 0], [0, -1/2]])
Q_vertical = as_tensor([[-1/2, 0], [0, 1/2]])
I = as_matrix([[1,0],[0,1]])
mat = grad(grad(u)) + q**2 * (Q+I/2) * u
E = (
+ a/2 * u**2 * dx
+ b/3 * u**3 * dx
+ c/4 * u**4 * dx
+ B * inner(mat, mat) * dx
+ K/2 * inner(grad(Q), grad(Q)) * dx
- l * tr(Q*Q) * dx
+ l * dot(tr(Q*Q), tr(Q*Q)) * dx
+ W/2 * inner(Q-Q_bottom, Q-Q_bottom) * ds(3)
+ W/2 * inner(Q-Q_vertical, Q-Q_vertical) * ds(4)
)
return E
def lagrangian(self, z, params):
E = self.energy(z, params)
return E
def residual(self, z, params, w):
L = self.lagrangian(z, params)
h = avg(CellDiameter(z.function_space().mesh()))
s = FacetNormal(z.function_space().mesh())
(u,_) = split(z)
(v,_) = split(w)
F = derivative(L, z, w) + h**(-3)*inner(jump(grad(u),s), jump(grad(v),s))*dS
return F
def boundary_conditions(self, Z, params):
return []
def initial_guess(self, Z, params, n):
# Tim's form for TFCD
(x, y) = SpatialCoordinate(Z.mesh())
R = Constant(1.5)
denomy = sqrt(R**2+x**2-2*R*sqrt(x**2)+y**2)+Constant(1e-10)
denomx = (sqrt(x**2)+Constant(1e-10)) * denomy
nx = x*(sqrt(x**2)-R)/denomy
ny = y/denomy
q0 = conditional(x**2>R**2, as_vector([-1/2, 0]),
as_vector([nx**2-1/2, nx*ny]))
lu = {"ksp_type": "preonly",
"pc_type": "lu",
"mat_type": "aij",
"pc_factor_mat_solver_type": "mumps",
"mat_mumps_icntl_14": 200,}
z = Function(Z)
z.sub(0).project(Constant(1.0), solver_parameters=lu)
z.sub(1).interpolate(q0)
return z
def solver_parameters(self, params, task, **kwargs):
params = {
"snes_max_it": 1000,
"snes_atol": 1.0e-8,
"snes_rtol": 1.0e-8,
"snes_monitor": None,
"snes_linesearch_type": "l2",
"snes_linesearch_monitor": None,
"snes_linesearch_maxstep": 1.0,
"snes_linesearch_damping": 1,
"snes_converged_reason": None,
"mat_type": "aij",
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"mat_mumps_icntl_14": 200,
"mat_mumps_icntl_24": 1,
"mat_mumps_icntl_13": 1
}
return params
def save_pvd(self, z, pvd, params):
mesh = z.function_space().mesh()
r = params[2]
(u, d) = z.split()
u.rename("Density")
#visualize the director
d0 = d[0]
d1 = d[1]
Q = interpolate(as_tensor([[d0, d1], [d1, -d0]]), TensorFunctionSpace(mesh, "CG", 1))
eigs, eigv = np.linalg.eigh(np.array(Q.vector()))
s = Function(FunctionSpace(mesh, "CG", 1)).interpolate(2*sqrt(dot(d,d)))
s.rename("order-parameter")
s_eig = Function(FunctionSpace(mesh, "CG", 1))
s_eig.vector()[:] = 2*eigs[:,1]
s_eig.rename("order-parameter-via-eig")
n = Function(VectorFunctionSpace(mesh, "CG", 1))
n.vector()[:,:] = eigv[:,:,1]
n.rename("Director")
pvd.write(u, n, s, s_eig)
def solver(self, problem, params, solver_params, prefix="", **kwargs):
Z = problem.u.function_space()
Xscat = self.Xscat
sign = self.sign
scatter = self.scatter
def post_function_callback(X, F):
# Basic idea: set residual to be value of slave dof + sign * value of master dof
self.scatter(X, Xscat, addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
with sign.dat.vec_ro as signv:
for slave_dof in self.my_slave:
#print(f"{Z.mesh().comm.rank}: setting value of {slave_dof} to be {X.getValue(slave_dof)} + {signv.getValue(slave_dof)} * {Xscat.getValue(slave_dof)} = {X.getValue(slave_dof) + signv.getValue(slave_dof) * Xscat.getValue(slave_dof)}", flush=True)
F.setValue(slave_dof, X.getValue(slave_dof) + signv.getValue(slave_dof) * Xscat.getValue(slave_dof))
F.assemble()
def post_jacobian_callback(X, J):
# Trust me. I know what I'm doing!
J.setOption(PETSc.Mat.Option.NEW_NONZERO_LOCATION_ERR, False)
J.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, False)
J.setOption(PETSc.Mat.Option.NEW_NONZERO_LOCATIONS, True)
J.setOption(PETSc.Mat.Option.UNUSED_NONZERO_LOCATION_ERR, False)
# First, zero all rows associated with slave dofs
J.zeroRows(self.my_slave, diag=1)
# Now set the off-diagonal entries
with sign.dat.vec_ro as signv:
for (slave_dof, master_dof) in zip(self.my_slave, self.my_master):
row = slave_dof
col = master_dof
val = signv.getValue(slave_dof)
J.setValues([row], [col], [val])
J.assemble()
solv = NonlinearVariationalSolver(problem, options_prefix=prefix,
solver_parameters=solver_params,
post_function_callback=post_function_callback,
post_jacobian_callback=post_jacobian_callback,
**kwargs)
return solv
if __name__ == "__main__":
problem = SmecticProblem()
comm = MPI.COMM_WORLD
mesh = problem.mesh(comm)
Z = problem.function_space(mesh)
params = [Constant(30), Constant(10), Constant(1)] # q, W, ratio
z = problem.initial_guess(Z, params, 0)
v = TestFunction(Z)
w = TrialFunction(Z)
F = problem.residual(z, params, v)
J = derivative(F, z, w)
bcs = problem.boundary_conditions(Z, params)
sp = problem.solver_parameters(params, None)
nproblem = NonlinearVariationalProblem(F, z, bcs=bcs, J=J)
solver = problem.solver(nproblem, params, sp)
solver.solve()
pvd = File("oilystreaks.pvd")
problem.save_pvd(z, pvd, params)I have given this quite a bit of thought and I think the most general case would be to have the option of periodic and aperiodic function spaces. I foresee problems with periodic topology in the event one would want to use a mixed function space where one of the spaces is aperiodic. I can't think of an example and it does seem like an edge case (I'll be here all week). I wouldn't know where to get started with this.
Consequently, I have decided to try to add support for periodic topologies. I have taken most of the code in firedrake-periodicity and have managed to get it to work in parallel. I am currently stuck on setting the coordinates on the new mesh. My idea is to mark the cells that contain vertices that have been remapped. I would like to construct a Coordinateless DG space by iterating over all cells and then finalize the mesh with this function. I am a little stuck with this step and further guidance would be appreciated. Adding MG support will be next.
` import firedrake.function as function import firedrake.functionspace as functionspace
dim=geometric_dim)
#coordinates_data = dmcommon.reordered_coords(plex, coordinates_fs.dm.getDefaultSection(), (self.num_vertices(), geometric_dim))
#coordinates = function.CoordinatelessFunction(coordinates_fs, val=coordinates_data, name="Coordinates")
vdg = functionspace.VectorFunctionSpace(self.topology, "DG", 1)
# HERE is where I would like to obtain the coord_data
coordinates_d = function.CoordinatelessFunction(vdg, val = coordinates_data)`I decided to try and give this a go from a different angle.
Let PETSc import the periodic topology and then use the multi-valued coordinate data to create CoordinatelessFunction that lives in a DG space.
Fortunately it seems to be working in serial but it fails in parallel with the following error: error.txt
Not sure how to debug this. Any advice is greatly appreciated.
wence-
commented
3 years ago This error is complaining that the coordinate array you are passing in doesn't have the right shape (number of entries). If things are working correctly in serial but not parallel, I suspect that you're not providing the halo data. Can you show your implementation?
medinaeder
commented
3 years ago I have added this to mesh.py In short I load in the periodic plex and then try to loop over the cells and just assume that the ordering works as a first pass. I am using what I would assume to be the standard distribution parameters. Perhaps this need to be updated for a DG function space?
`def callback(self):
del self._callback
# Finish the initialisation of mesh topology
self.topology.init()
coordinates_fs = functionspace.VectorFunctionSpace(self.topology, "Lagrange", 1,
dim=geometric_dim)
if True:
coordinates_fs = functionspace.VectorFunctionSpace(self.topology, "Discontinuous Lagrange", 1,
dim=geometric_dim)
dmp = self._topology_dm
dm_cells = dmp.getHeightStratum(0)
print(dm_cells)
dm_vertices = dmp.getDepthStratum(0)
dm_coord_sec = dmp.getCoordinateSection()
dm_section = coordinates_fs.dm.getDefaultSection()
#dm_coordinates = plex.getCoordinatesLocal().array.reshape((-1,geometric_dim))
coordinates_data = np.empty(coordinates_fs.dim()).reshape((-1,geometric_dim))
#import matplotlib.pyplot as plt
for p in range(*dm_cells):
x = dmp.getVecClosure(dm_coord_sec, dmp.getCoordinatesLocal(),p).reshape((-1,2))
#plt.scatter(x[:3,0], x[:3,1])
# FIXME:
spacedim = 2 # Not asssuming embeded surfaces
for i in range(spacedim+1):
for j in range(spacedim):
coordinates_data[3*p+i,j] = x[i,j]
coordinates = function.CoordinatelessFunction(coordinates_fs,
val=coordinates_data,
name="Coordinates")
#import IPython; IPython.embed()
else:
coordinates_data = dmcommon.reordered_coords(plex, coordinates_fs.dm.getDefaultSection(),
(self.num_vertices(), geometric_dim))
coordinates = function.CoordinatelessFunction(coordinates_fs,
val=coordinates_data,
name="Coordinates")
self.__init__(coordinates)
mesh._callback = callback
return mesh`
coordinates_data = np.empty(coordinates_fs.dim()).reshape((-1,geometric_dim))
This line is the first problem I see. coordinates_fs.dim() is the global size of the functionspace. You can instead use.
coordinates_data = np.empty((coordinates_fs.node_set.total_size, geometric_dim))I think.
medinaeder
commented
3 years ago Hi everyone!
Some progress has been made. The code seems to be working in serial for both 2D and 3D geometries!
On each processor I loop through the cell_closure list, extract the multivalued coordinates and create two lists (if multivalued). I then pull the local node value and check to see if it is the the list that contains the periodic coordinates. I am assume that the ordering in the cell_closure is the ordering of the nodes on each element.
I forgot to pose the question:
Is this assumption valid? Is there an offset that i am not taking into account?
def Mesh(... , periodic=True):
...
def callback(self):
"""Finish initialisation."""
del self._callback
# Finish the initialisation of mesh topology
self.topology.init()
if self._periodic:
print("Attempting Periodic Conversion")
coordinates_fs = functionspace.VectorFunctionSpace(self.topology,
"Discontinuous Lagrange", 1,
dim=geometric_dim)
# So every proc has a set of mesh entities that have already been partitioned
fs_section = coordinates_fs.dm.getDefaultSection()
coordinates_data = np.empty((coordinates_fs.node_set.total_size,geometric_dim))
dm_cells = plex.getHeightStratum(0)
dm_coord_sec = plex.getCoordinateSection()
closure = self._topology.cell_closure
#with np.printoptions(threshold=np.inf):
# print(closure)
simplex_dim = geometric_dim+1
# Could be sped up only looping through cells became periodic
for row,cell in enumerate(closure[:,-1]):
ys = plex.getVecClosure(dm_coord_sec,
plex.getCoordinatesLocal(),cell).reshape((-1,geometric_dim))
y = ys[:simplex_dim,:] # Keep coords generated prior to periodization
y2 = ys[simplex_dim:2*simplex_dim,:] # Pull multivalued
x = np.empty_like(y)
for node in range(simplex_dim):
dm_node = closure[row,node]
x_temp = plex.getVecClosure(dm_coord_sec, plex.getCoordinatesLocal(),
dm_node).reshape((-1, geometric_dim))
# The node should be in one of the two lists.
missing_node = x_temp[0].tolist() not in y.tolist()
if missing_node:
loc = np.argmax(np.sum(x_temp == y2,axis=1)) # Find the closests point
x[node,:] = y[loc,:]
else:
x[node,:] = x_temp
# What could be causing the scramble in parallel?
for node in range(simplex_dim):
for j in range(geometric_dim):
# Need to get the correct ordering
coordinates_data[simplex_dim*cell+node,j] = x[node,j]
coordinates = function.CoordinatelessFunction(coordinates_fs,
val=coordinates_data,
name="Coordinates")
...
wence-
commented
3 years ago # What could be causing the scramble in parallel?
for node in range(simplex_dim):
for j in range(geometric_dim):
# Need to get the correct ordering
coordinates_data[simplex_dim*cell+node,j] = x[node,j]Perhaps
firedrake_coord_sec = coordinates_fs.dm.getLocalSection()
offset = firedrake_coord_sec.getOffset(cell)
dof = firedrake_coord_sec.getDof(cell)
assert dof == simplex_dim
coordinates_data[offset:offset+dof, :] = x[:]?
wence-
commented
3 years ago The other thing I can think is that we are assuming as well that the vertices are ordered in the same way in firedrake as plex (not sure this is true).
One final bit, you need to make a FiniteElement for this field with variant="equispaced" otherwise I think the DG nodes are not physically located at the vertices (which I suspect the mesh generator does)
wence-
commented
3 years ago Can you show the bad output VTU files on the two separate processes, that might give us a hint as well.
wence-
commented
3 years ago # What could be causing the scramble in parallel? for node in range(simplex_dim): for j in range(geometric_dim): # Need to get the correct ordering coordinates_data[simplex_dim*cell+node,j] = x[node,j]Perhaps
firedrake_coord_sec = coordinates_fs.dm.getLocalSection() offset = firedrake_coord_sec.getOffset(cell) dof = firedrake_coord_sec.getDof(cell) assert dof == simplex_dim coordinates_data[offset:offset+dof, :] = x[:]
This is slightly more efficient, but results in, I think, the same numbering, so it's probably not that.
pefarrell
commented
3 years ago Does @medinaeder 's code live on a branch somewhere? I looked but couldn't find it. It would be nice to try it (in serial).
medinaeder
commented
3 years ago It is on a local branch. Can I get push access to firedrake?
dham
commented
3 years ago It is on a local branch. Can I get push access to firedrake?
You should have an invite now.
medinaeder
commented
3 years ago Awesome thanks!
The branch is called em/periodic
gonzalogddiego
commented
3 years ago Hello everyone,
I am interested in using periodic meshes from gmsh and I have just tried the new branch from @medinaeder . However, if I am on the branch and I simply run
from firedrake import *
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 1)I get the following error from UFL:
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "<decorator-gen-24>", line 2, in FunctionSpace
File "/home/gonzalo/Documents/firedrake/src/PyOP2/pyop2/profiling.py", line 60, in wrapper
return f(*args, **kwargs)
File "/home/gonzalo/Documents/firedrake/src/firedrake/firedrake/functionspace.py", line 117, in FunctionSpace
element = make_scalar_element(mesh, family, degree, vfamily, vdegree)
File "/home/gonzalo/Documents/firedrake/src/firedrake/firedrake/functionspace.py", line 40, in make_scalar_element
mesh.init()
File "/home/gonzalo/Documents/firedrake/src/firedrake/firedrake/mesh.py", line 1342, in init
self._callback(self)
File "/home/gonzalo/Documents/firedrake/src/firedrake/firedrake/mesh.py", line 1813, in callback
dim=geometric_dim)
File "/home/gonzalo/Documents/firedrake/src/firedrake/firedrake/functionspace.py", line 180, in VectorFunctionSpace
element = ufl.VectorElement(sub_element, dim=dim, variant=variant)
TypeError: __init__() got an unexpected keyword argument 'variantDo I need a particular branch of UFL to get this working? I updated firedrake just before fetching this branch.
Thanks!
medinaeder
commented
3 years ago I did change my UFL locally to accept the variant keyword but I don't know if it is truly needed.
Also note that the branch is not stable. I made a mistake when I said it worked in serial. Apparently it only worked in serial when solving poissons equation on the periodic mesh that I had created. You can try this out for yourself. It is under tests/periodic/. When I tried to solve a linear elasticity problem it seems like it is treating it as a discontinuous problem.
I think I know the fix. My apologies for the inconvenience.
dthillaithevan
commented
1 year ago Hey, just wanted to see if there has been progress on this? I am hoping to solve periodic, non-linear elasticity problems using imported periodic meshes (generated in GMSH).
martintruffer
commented
2 months ago Same question here.
Dear Firedrake developers,
I am currently working on some simulations using periodic boundary conditions. I thought it would be a good idea to summarize all information that I could find up to now and trigger a discussion about the implementation:
PeriodicSquareMeshor thePeriodicBoxMesh(see #1798 ).1246 suggests that petsc is able to ready periodic meshes, but one needs to create a proper discontinuous coordinate field.
dolfinone just uses a non-periodic mesh and then defines a master-slave map, which eliminates matching boundary nodes. This map is passed during the creation of theFunctionSpaceFrom my point of view the
dolfinversion has some advantages:dolfincode is available and can perhaps be adapted to firedrakeI write that the
dolfincode can be adapted, but unfortunately I was not able to do it by myself. At least not without some additional inputs. As far as I could see the original code was commited by Garth N. Wells. Maybe he could tell how much effort it would be to adapt the code (and how to do it :-)What do you think?