Loading of periodic meshes generated by gmsh

firedrake/mesh.py

@@ -3190,6 +3190,76 @@ def make_mesh_from_coordinates(coordinates, name, tolerance=0.5):
     return mesh
 
 
+def _fully_localize_coordinates(dm: PETSc.DM) -> None:
+    """Expand partially localized coordinates to cover all cells.
+
+    PETSc's ``DMLocalizeCoordinates`` only creates cell-local (DG)
+    coordinates for cells that straddle the periodic boundary.  The
+    Firedrake DG coordinate path in ``reordered_coords`` expects every
+    cell to have an entry in the cell coordinate vector.  This helper
+    fills in the missing entries using the CG vertex coordinates
+    obtained via ``vecGetClosure``.
+
+    Parameters
+    ----------
+    dm : PETSc.DM
+        The DM whose cell coordinates should be fully localized.
+    """
+    if not dm.getCoordinatesLocalized():
+        return
+
+    gdim = dm.getCoordinateDim()
+    cStart, cEnd = dm.getHeightStratum(0)
+    cell_sec = dm.getCellCoordinateSection()
+
+    dofs_per_cell = None
+    needs_expansion = False
+    for c in range(cStart, cEnd):
+        dof = cell_sec.getDof(c)
+        if dof > 0:
+            dofs_per_cell = dof
+        else:
+            needs_expansion = True
+
+    if not needs_expansion or dofs_per_cell is None:
+        return
+
+    coord_sec = dm.getCoordinateSection()
+    coord_vec = dm.getCoordinatesLocal()
+    old_cell_vec = dm.getCellCoordinatesLocal()
+
+    new_sec = PETSc.Section().create(comm=PETSc.COMM_SELF)
+    new_sec.setNumFields(1)
+    new_sec.setFieldComponents(0, gdim)
+    new_sec.setChart(cStart, cEnd)
+    for c in range(cStart, cEnd):
+        new_sec.setDof(c, dofs_per_cell)
+        new_sec.setFieldDof(c, 0, dofs_per_cell)
+    new_sec.setUp()
+
+    new_vec = PETSc.Vec().create(comm=PETSc.COMM_SELF)
+    new_vec.setSizes((new_sec.getStorageSize(), PETSc.DETERMINE), gdim)
+    new_vec.setType(coord_vec.getType())
+
+    arr = new_vec.array
+    for c in range(cStart, cEnd):
+        new_offset = new_sec.getOffset(c)
+        old_dof = cell_sec.getDof(c)
+        if old_dof > 0:
+            old_offset = cell_sec.getOffset(c)
+            arr[new_offset:new_offset + dofs_per_cell] = \
+                old_cell_vec.array[old_offset:old_offset + old_dof]
+        else:
+            closure = dm.vecGetClosure(coord_sec, coord_vec, c)
+            arr[new_offset:new_offset + dofs_per_cell] = closure[:dofs_per_cell]
+
+    coord_dm = dm.getCoordinateDM()
+    new_coord_dm = coord_dm.clone()
+    dm.setCellCoordinateDM(new_coord_dm)
+    dm.setCellCoordinateSection(gdim, new_sec)
+    dm.setCellCoordinatesLocal(new_vec)
+
+
 def make_mesh_from_mesh_topology(topology, name, tolerance=0.5):
     """Make mesh from topology.
 
@@ -3212,6 +3282,7 @@ def make_mesh_from_mesh_topology(topology, name, tolerance=0.5):
     # TODO: meshfile might indicates higher-order coordinate element
     cell = topology.ufl_cell()
     geometric_dim = topology.topology_dm.getCoordinateDim()
+    _fully_localize_coordinates(topology.topology_dm)
     if not topology.topology_dm.getCoordinatesLocalized():
         element = finat.ufl.VectorElement("Lagrange", cell, 1, dim=geometric_dim)
     else:

tests/periodic/geom/p2d.geo

@@ -0,0 +1,16 @@
+SetFactory("OpenCASCADE");
+//+
+Rectangle(1) = {0, 0, 0, 0.6, 0.5, 0};
+MeshSize {:} = 0.1;
+MeshSize {:} = 0.1;
+
+// S2: Right Side
+// S4: Left Side
+
+Periodic Curve {2} = {4} Translate {0.6,0,0};
+
+Physical Surface(1) = {1};
+Physical Curve(1) = {1};
+Physical Curve(2) = {2};
+Physical Curve(3) = {3};
+Physical Curve(4) = {4};

tests/periodic/geom/p2d.msh

@@ -0,0 +1,175 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+52
+1 0 0 0
+2 0.6 0 0
+3 0.6 0.5 0
+4 0 0.5 0
+5 0.1 0 0
+6 0.2000000000000002 0 0
+7 0.3000000000000003 0 0
+8 0.4000000000000006 0 0
+9 0.5000000000000007 0 0
+10 0.6 0.09999999999999998 0
+11 0.6 0.2 0
+12 0.6 0.3 0
+13 0.6 0.4 0
+14 0.4999999999999999 0.5 0
+15 0.3999999999999997 0.5 0
+16 0.2999999999999997 0.5 0
+17 0.1999999999999993 0.5 0
+18 0.09999999999999931 0.5 0
+19 0 0.4 0
+20 0 0.3 0
+21 0 0.2 0
+22 0 0.09999999999999998 0
+23 0.5186224219100395 0.245166047252571 0
+24 0.08137757808996043 0.2451660472525709 0
+25 0.2543372475932793 0.08224840025510172 0
+26 0.3499999999999996 0.413397459621556 0
+27 0.4430574538662087 0.09595102731164898 0
+28 0.150423736928087 0.415322795494417 0
+29 0.2500706228213473 0.4137183489336996 0
+30 0.3000117704702239 0.3268484007951361 0
+31 0.4018525930917393 0.3349783789207162 0
+32 0.3578740082166879 0.2419001089813815 0
+33 0.2435188692474362 0.2409750047583969 0
+34 0.4497777049626761 0.41512773648052 0
+35 0.3452285359833018 0.08308079139007907 0
+36 0.1575621529328313 0.09532900729402857 0
+37 0.1971535352311468 0.3354831874338384 0
+38 0.08080880024257878 0.1530389498800314 0
+39 0.5059889644735611 0.3418579238694238 0
+40 0.5204794770710884 0.1523966756079144 0
+41 0.0938978312712459 0.3419579049447122 0
+42 0.4431474979352846 0.1938733509157883 0
+43 0.1574026289941869 0.1937225425877068 0
+44 0.2999999999999989 0.1633974596215566 0
+45 0.5267949192431121 0.4267949192431121 0
+46 0.07320508075688786 0.07320508075688781 0
+47 0.5267949192431126 0.0732050807568875 0
+48 0.07320508075688736 0.4267949192431125 0
+49 0.3759586495219419 0.156251476709271 0
+50 0.2251229998750313 0.1564314424076824 0
+51 0.4454970971254625 0.2715551619879761 0
+52 0.1551248878926393 0.2709383851890471 0
+$EndNodes
+$Elements
+102
+1 1 2 1 1 1 5
+2 1 2 1 1 5 6
+3 1 2 1 1 6 7
+4 1 2 1 1 7 8
+5 1 2 1 1 8 9
+6 1 2 1 1 9 2
+7 1 2 2 2 2 10
+8 1 2 2 2 10 11
+9 1 2 2 2 11 12
+10 1 2 2 2 12 13
+11 1 2 2 2 13 3
+12 1 2 3 3 3 14
+13 1 2 3 3 14 15
+14 1 2 3 3 15 16
+15 1 2 3 3 16 17
+16 1 2 3 3 17 18
+17 1 2 3 3 18 4
+18 1 2 4 4 4 19
+19 1 2 4 4 19 20
+20 1 2 4 4 20 21
+21 1 2 4 4 21 22
+22 1 2 4 4 22 1
+23 2 2 1 1 31 39 34
+24 2 2 1 1 28 41 37
+25 2 2 1 1 40 42 27
+26 2 2 1 1 36 43 38
+27 2 2 1 1 31 51 39
+28 2 2 1 1 41 52 37
+29 2 2 1 1 42 49 27
+30 2 2 1 1 36 50 43
+31 2 2 1 1 15 34 14
+32 2 2 1 1 29 30 26
+33 2 2 1 1 7 35 25
+34 2 2 1 1 16 29 26
+35 2 2 1 1 30 37 33
+36 2 2 1 1 9 27 8
+37 2 2 1 1 29 37 30
+38 2 2 1 1 6 36 5
+39 2 2 1 1 26 34 15
+40 2 2 1 1 8 35 7
+41 2 2 1 1 17 29 16
+42 2 2 1 1 18 28 17
+43 2 2 1 1 30 31 26
+44 2 2 1 1 30 32 31
+45 2 2 1 1 25 36 6
+46 2 2 1 1 16 26 15
+47 2 2 1 1 7 25 6
+48 2 2 1 1 21 24 20
+49 2 2 1 1 31 34 26
+50 2 2 1 1 28 37 29
+51 2 2 1 1 30 33 32
+52 2 2 1 1 12 23 11
+53 2 2 1 1 27 35 8
+54 2 2 1 1 28 29 17
+55 2 2 1 1 22 38 21
+56 2 2 1 1 13 39 12
+57 2 2 1 1 11 40 10
+58 2 2 1 1 20 41 19
+59 2 2 1 1 21 38 24
+60 2 2 1 1 23 40 11
+61 2 2 1 1 24 41 20
+62 2 2 1 1 12 39 23
+63 2 2 1 1 1 46 22
+64 2 2 1 1 2 47 9
+65 2 2 1 1 3 45 13
+66 2 2 1 1 5 46 1
+67 2 2 1 1 14 45 3
+68 2 2 1 1 10 47 2
+69 2 2 1 1 19 48 4
+70 2 2 1 1 4 48 18
+71 2 2 1 1 27 49 35
+72 2 2 1 1 25 50 36
+73 2 2 1 1 35 44 25
+74 2 2 1 1 33 44 32
+75 2 2 1 1 23 42 40
+76 2 2 1 1 38 43 24
+77 2 2 1 1 32 51 31
+78 2 2 1 1 37 52 33
+79 2 2 1 1 39 51 23
+80 2 2 1 1 24 52 41
+81 2 2 1 1 42 51 32
+82 2 2 1 1 33 52 43
+83 2 2 1 1 41 48 19
+84 2 2 1 1 40 47 10
+85 2 2 1 1 22 46 38
+86 2 2 1 1 13 45 39
+87 2 2 1 1 34 45 14
+88 2 2 1 1 36 46 5
+89 2 2 1 1 9 47 27
+90 2 2 1 1 18 48 28
+91 2 2 1 1 27 47 40
+92 2 2 1 1 28 48 41
+93 2 2 1 1 38 46 36
+94 2 2 1 1 39 45 34
+95 2 2 1 1 43 50 33
+96 2 2 1 1 32 49 42
+97 2 2 1 1 33 50 44
+98 2 2 1 1 44 49 32
+99 2 2 1 1 23 51 42
+100 2 2 1 1 43 52 24
+101 2 2 1 1 35 49 44
+102 2 2 1 1 44 50 25
+$EndElements
+$Periodic
+1
+1 2 4
+Affine 1 0 0 0.6 0 1 0 0 0 0 1 0 0 0 0 1
+6
+2 1
+3 4
+10 22
+11 21
+12 20
+13 19
+$EndPeriodic

tests/periodic/geom/p2d_xy.geo

@@ -0,0 +1,16 @@
+SetFactory("OpenCASCADE");
+//+
+Rectangle(1) = {0, 0, 0, 0.6, 0.5, 0};
+MeshSize {:} = 0.05;
+
+// Curve 1: bottom (y=0), Curve 2: right (x=0.6)
+// Curve 3: top (y=0.5), Curve 4: left (x=0)
+
+Periodic Curve {2} = {4} Translate {0.6, 0, 0};
+Periodic Curve {3} = {1} Translate {0, 0.5, 0};
+
+Physical Surface(1) = {1};
+Physical Curve(1) = {1};
+Physical Curve(2) = {2};
+Physical Curve(3) = {3};
+Physical Curve(4) = {4};