Skip to content

Commit

Permalink
Extend to 2D-0D and use updated UFL nomenclature
Browse files Browse the repository at this point in the history
  • Loading branch information
jorgensd committed Jan 14, 2025
1 parent 5fa9d30 commit 7a9c0cf
Show file tree
Hide file tree
Showing 12 changed files with 265 additions and 131 deletions.
2 changes: 1 addition & 1 deletion ffcx/codegeneration/C/form.py
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ def generator(ir: FormIR, options):
integral_ids = []
integral_offsets = [0]
# Note: the order of this list is defined by the enum ufcx_integral_type in ufcx.h
for itg_type in ("cell", "exterior_facet", "interior_facet", "edge"):
for itg_type in ("cell", "exterior_facet", "interior_facet", "ridge"):
unsorted_integrals = []
unsorted_ids = []
for name, id in zip(ir.integral_names[itg_type], ir.subdomain_ids[itg_type]):
Expand Down
67 changes: 34 additions & 33 deletions ffcx/codegeneration/access.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,14 +40,14 @@ def __init__(self, entity_type: str, integral_type: str, symbols, options):
ufl.geometry.CellCoordinate: self.cell_coordinate,
ufl.geometry.FacetCoordinate: self.facet_coordinate,
ufl.geometry.CellVertices: self.cell_vertices,
ufl.geometry.FacetEdgeVectors: self.facet_edge_vectors,
ufl.geometry.CellEdgeVectors: self.cell_edge_vectors,
ufl.geometry.FacetEdgeVectors: self.facet_ridge_vectors,
ufl.geometry.CellEdgeVectors: self.cell_ridge_vectors,
ufl.geometry.CellFacetJacobian: self.cell_facet_jacobian,
ufl.geometry.CellEdgeJacobian: self.cell_edge_jacobian,
ufl.geometry.CellRidgeJacobian: self.cell_ridge_jacobian,
ufl.geometry.ReferenceCellVolume: self.reference_cell_volume,
ufl.geometry.ReferenceFacetVolume: self.reference_facet_volume,
ufl.geometry.ReferenceCellEdgeVectors: self.reference_cell_edge_vectors,
ufl.geometry.ReferenceFacetEdgeVectors: self.reference_facet_edge_vectors,
ufl.geometry.ReferenceCellEdgeVectors: self.reference_cell_ridge_vectors,
ufl.geometry.ReferenceFacetEdgeVectors: self.reference_facet_ridge_vectors,
ufl.geometry.ReferenceNormal: self.reference_normal,
ufl.geometry.CellOrientation: self._pass,
ufl.geometry.FacetOrientation: self.facet_orientation,
Expand All @@ -71,7 +71,6 @@ def get(
if isinstance(e, k):
handler = self.call_lookup[k]
break

if handler:
return handler(mt, tabledata, quadrature_rule) # type: ignore
else:
Expand Down Expand Up @@ -247,38 +246,40 @@ def cell_facet_jacobian(self, mt, tabledata, num_points):
else:
raise RuntimeError(f"Unhandled cell types {cellname}.")

def cell_edge_jacobian(self, mt, tabledata, num_points):
"""Access a cell edge jacobian."""
def cell_ridge_jacobian(self, mt, tabledata, num_points):
"""Access a cell ridge jacobian."""
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname in ("tetrahedron", "prism", "hexahedron"):
table = L.Symbol(f"{cellname}_cell_edge_jacobian", dtype=L.DataType.REAL)
edge = self.symbols.entity("edge", mt.restriction)
return table[edge][mt.component[0]][mt.component[1]]
table = L.Symbol(f"{cellname}_cell_ridge_jacobian", dtype=L.DataType.REAL)
ridge = self.symbols.entity("ridge", mt.restriction)
return table[ridge][mt.component[0]][mt.component[1]]
elif cellname in ["triangle", "quadrilateral"]:
raise RuntimeError("The reference facet jacobian doesn't make sense for 2D cells.")
else:
raise RuntimeError(f"Unhandled cell types {cellname}.")

def reference_cell_edge_vectors(self, mt, tabledata, num_points):
"""Access a reference cell edge vector."""
def reference_cell_ridge_vectors(self, mt, tabledata, num_points):
"""Access a reference cell ridge vector."""
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname in ("triangle", "tetrahedron", "quadrilateral", "hexahedron"):
table = L.Symbol(f"{cellname}_reference_cell_edge_vectors", dtype=L.DataType.REAL)
table = L.Symbol(f"{cellname}_reference_cell_ridge_vectors", dtype=L.DataType.REAL)
return table[mt.component[0]][mt.component[1]]
elif cellname == "interval":
raise RuntimeError(
"The reference cell edge vectors doesn't make sense for interval cell."
"The reference cell ridge vectors doesn't make sense for interval cell."
)
else:
raise RuntimeError(f"Unhandled cell types {cellname}.")

def reference_facet_edge_vectors(self, mt, tabledata, num_points):
"""Access a reference facet edge vector."""
def reference_facet_ridge_vectors(self, mt, tabledata, num_points):
"""Access a reference facet ridge vector."""
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname in ("tetrahedron", "hexahedron"):
table = L.Symbol(f"{cellname}_reference_facet_edge_vectors", dtype=L.DataType.REAL)
table = L.Symbol(f"{cellname}_reference_facet_ridge_vectors", dtype=L.DataType.REAL)
return table[mt.component[0]][mt.component[1]]
elif cellname in ("interval", "triangle", "quadrilateral"):
raise RuntimeError(
"The reference cell facet edge vectors doesn't make sense for interval "
"The reference cell facet ridge vectors doesn't make sense for interval "
"or triangle cell."
)
else:
Expand Down Expand Up @@ -318,8 +319,8 @@ def cell_vertices(self, mt, tabledata, num_points):
expr = self.symbols.domain_dof_access(dof, component, gdim, num_scalar_dofs, mt.restriction)
return expr

def cell_edge_vectors(self, mt, tabledata, num_points):
"""Access a cell edge vector."""
def cell_ridge_vectors(self, mt, tabledata, num_points):
"""Access a cell ridge vector."""
# Get properties of domain
domain = ufl.domain.extract_unique_domain(mt.terminal)
cellname = domain.ufl_cell().cellname()
Expand All @@ -330,7 +331,7 @@ def cell_edge_vectors(self, mt, tabledata, num_points):
pass
elif cellname == "interval":
raise RuntimeError(
"The physical cell edge vectors doesn't make sense for interval cell."
"The physical cell ridge vectors doesn't make sense for interval cell."
)
else:
raise RuntimeError(f"Unhandled cell types {cellname}.")
Expand All @@ -345,9 +346,9 @@ def cell_edge_vectors(self, mt, tabledata, num_points):
vertex_scalar_dofs = scalar_element.entity_dofs[0]
num_scalar_dofs = scalar_element.dim

# Get edge vertices
edge = mt.component[0]
vertex0, vertex1 = scalar_element.reference_topology[1][edge]
# Get ridge vertices
ridge = mt.component[0]
vertex0, vertex1 = scalar_element.reference_topology[1][ridge]

# Get dofs and component
(dof0,) = vertex_scalar_dofs[vertex0]
Expand All @@ -358,8 +359,8 @@ def cell_edge_vectors(self, mt, tabledata, num_points):
dof0, component, gdim, num_scalar_dofs, mt.restriction
) - self.symbols.domain_dof_access(dof1, component, gdim, num_scalar_dofs, mt.restriction)

def facet_edge_vectors(self, mt, tabledata, num_points):
"""Access a facet edge vector."""
def facet_ridge_vectors(self, mt, tabledata, num_points):
"""Access a facet ridge vector."""
# Get properties of domain
domain = ufl.domain.extract_unique_domain(mt.terminal)
cellname = domain.ufl_cell().cellname()
Expand All @@ -370,7 +371,7 @@ def facet_edge_vectors(self, mt, tabledata, num_points):
pass
elif cellname in ("interval", "triangle", "quadrilateral"):
raise RuntimeError(
f"The physical facet edge vectors doesn't make sense for {cellname} cell."
f"The physical facet ridge vectors doesn't make sense for {cellname} cell."
)
else:
raise RuntimeError(f"Unhandled cell types {cellname}.")
Expand All @@ -385,12 +386,12 @@ def facet_edge_vectors(self, mt, tabledata, num_points):
scalar_element = ufl_scalar_element
num_scalar_dofs = scalar_element.dim

# Get edge vertices
# Get ridge vertices
facet = self.symbols.entity("facet", mt.restriction)
facet_edge = mt.component[0]
facet_edge_vertices = L.Symbol(f"{cellname}_facet_edge_vertices", dtype=L.DataType.INT)
vertex0 = facet_edge_vertices[facet][facet_edge][0]
vertex1 = facet_edge_vertices[facet][facet_edge][1]
facet_ridge = mt.component[0]
facet_ridge_vertices = L.Symbol(f"{cellname}_facet_ridge_vertices", dtype=L.DataType.INT)
vertex0 = facet_ridge_vertices[facet][facet_ridge][0]
vertex1 = facet_ridge_vertices[facet][facet_ridge][1]

# Get dofs and component
component = mt.component[1]
Expand Down
2 changes: 1 addition & 1 deletion ffcx/codegeneration/definitions.py
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ def __init__(self, entity_type: str, integral_type: str, access, options):
ufl.geometry.FacetEdgeVectors: self.pass_through,
ufl.geometry.CellEdgeVectors: self.pass_through,
ufl.geometry.CellFacetJacobian: self.pass_through,
ufl.geometry.CellEdgeJacobian: self.pass_through,
ufl.geometry.CellRidgeJacobian: self.pass_through,
ufl.geometry.ReferenceCellVolume: self.pass_through,
ufl.geometry.ReferenceFacetVolume: self.pass_through,
ufl.geometry.ReferenceCellEdgeVectors: self.pass_through,
Expand Down
6 changes: 3 additions & 3 deletions ffcx/codegeneration/expression_generator.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,12 +62,12 @@ def generate(self):
def generate_geometry_tables(self):
"""Generate static tables of geometry data."""
ufl_geometry = {
ufl.geometry.FacetEdgeVectors: "facet_edge_vectors",
ufl.geometry.FacetEdgeVectors: "facet_ridge_vectors",
ufl.geometry.CellFacetJacobian: "cell_facet_jacobian",
ufl.geometry.ReferenceCellVolume: "reference_cell_volume",
ufl.geometry.ReferenceFacetVolume: "reference_facet_volume",
ufl.geometry.ReferenceCellEdgeVectors: "reference_cell_edge_vectors",
ufl.geometry.ReferenceFacetEdgeVectors: "reference_facet_edge_vectors",
ufl.geometry.ReferenceCellEdgeVectors: "reference_cell_ridge_vectors",
ufl.geometry.ReferenceFacetEdgeVectors: "reference_facet_ridge_vectors",
ufl.geometry.ReferenceNormal: "reference_normals",
}

Expand Down
64 changes: 32 additions & 32 deletions ffcx/codegeneration/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,47 +13,47 @@

def write_table(tablename, cellname):
"""Write a table."""
if tablename == "facet_edge_vertices":
return facet_edge_vertices(tablename, cellname)
if tablename == "facet_ridge_vertices":
return facet_ridge_vertices(tablename, cellname)
if tablename == "cell_facet_jacobian":
return cell_facet_jacobian(tablename, cellname)
if tablename == "cell_edge_jacobian":
return cell_edge_jacobian(tablename, cellname)
if tablename == "cell_ridge_jacobian":
return cell_ridge_jacobian(tablename, cellname)
if tablename == "reference_cell_volume":
return reference_cell_volume(tablename, cellname)
if tablename == "reference_facet_volume":
return reference_facet_volume(tablename, cellname)
if tablename == "reference_cell_edge_vectors":
return reference_cell_edge_vectors(tablename, cellname)
if tablename == "reference_facet_edge_vectors":
return reference_facet_edge_vectors(tablename, cellname)
if tablename == "reference_cell_ridge_vectors":
return reference_cell_ridge_vectors(tablename, cellname)
if tablename == "reference_facet_ridge_vectors":
return reference_facet_ridge_vectors(tablename, cellname)
if tablename == "reference_normals":
return reference_normals(tablename, cellname)
if tablename == "facet_orientation":
return facet_orientation(tablename, cellname)
raise ValueError(f"Unknown geometry table name: {tablename}")


def facet_edge_vertices(tablename, cellname):
"""Write facet edge vertices."""
def facet_ridge_vertices(tablename, cellname):
"""Write facet ridge vertices."""
celltype = getattr(basix.CellType, cellname)
topology = basix.topology(celltype)
triangle_edges = basix.topology(basix.CellType.triangle)[1]
quadrilateral_edges = basix.topology(basix.CellType.quadrilateral)[1]
triangle_ridges = basix.topology(basix.CellType.triangle)[1]
quadrilateral_ridges = basix.topology(basix.CellType.quadrilateral)[1]

if len(topology) != 4:
raise ValueError("Can only get facet edges for 3D cells.")
raise ValueError("Can only get facet ridges for 3D cells.")

edge_vertices = []
ridge_vertices = []
for facet in topology[-2]:
if len(facet) == 3:
edge_vertices += [[[facet[i] for i in edge] for edge in triangle_edges]]
ridge_vertices += [[[facet[i] for i in ridge] for ridge in triangle_ridges]]
elif len(facet) == 4:
edge_vertices += [[[facet[i] for i in edge] for edge in quadrilateral_edges]]
ridge_vertices += [[[facet[i] for i in ridge] for ridge in quadrilateral_ridges]]
else:
raise ValueError("Only triangular and quadrilateral faces supported.")

out = np.array(edge_vertices, dtype=int)
out = np.array(ridge_vertices, dtype=int)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.INT)
return L.ArrayDecl(symbol, values=out, const=True)

Expand All @@ -66,7 +66,7 @@ def cell_facet_jacobian(tablename, cellname):
return L.ArrayDecl(symbol, values=out, const=True)


def cell_edge_jacobian(tablename, cellname):
def cell_ridge_jacobian(tablename, cellname):
"""Write a reference facet jacobian."""
celltype = getattr(basix.CellType, cellname)
out = basix.cell.edge_jacobians(celltype)
Expand All @@ -93,40 +93,40 @@ def reference_facet_volume(tablename, cellname):
return L.VariableDecl(symbol, volumes[0])


def reference_cell_edge_vectors(tablename, cellname):
"""Write reference edge vectors."""
def reference_cell_ridge_vectors(tablename, cellname):
"""Write reference ridge vectors."""
celltype = getattr(basix.CellType, cellname)
topology = basix.topology(celltype)
geometry = basix.geometry(celltype)
edge_vectors = [geometry[j] - geometry[i] for i, j in topology[1]]
out = np.array(edge_vectors)
ridge_vectors = [geometry[j] - geometry[i] for i, j in topology[1]]
out = np.array(ridge_vectors)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.ArrayDecl(symbol, values=out, const=True)


def reference_facet_edge_vectors(tablename, cellname):
"""Write facet reference edge vectors."""
def reference_facet_ridge_vectors(tablename, cellname):
"""Write facet reference ridge vectors."""
celltype = getattr(basix.CellType, cellname)
topology = basix.topology(celltype)
geometry = basix.geometry(celltype)
triangle_edges = basix.topology(basix.CellType.triangle)[1]
quadrilateral_edges = basix.topology(basix.CellType.quadrilateral)[1]
triangle_ridges = basix.topology(basix.CellType.triangle)[1]
quadrilateral_ridges = basix.topology(basix.CellType.quadrilateral)[1]

if len(topology) != 4:
raise ValueError("Can only get facet edges for 3D cells.")
raise ValueError("Can only get facet ridges for 3D cells.")

edge_vectors = []
ridge_vectors = []
for facet in topology[-2]:
if len(facet) == 3:
edge_vectors += [geometry[facet[j]] - geometry[facet[i]] for i, j in triangle_edges]
ridge_vectors += [geometry[facet[j]] - geometry[facet[i]] for i, j in triangle_ridges]
elif len(facet) == 4:
edge_vectors += [
geometry[facet[j]] - geometry[facet[i]] for i, j in quadrilateral_edges
ridge_vectors += [
geometry[facet[j]] - geometry[facet[i]] for i, j in quadrilateral_ridges
]
else:
raise ValueError("Only triangular and quadrilateral faces supported.")

out = np.array(edge_vectors)
out = np.array(ridge_vectors)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.ArrayDecl(symbol, values=out, const=True)

Expand Down
2 changes: 1 addition & 1 deletion ffcx/codegeneration/integral_generator.py
Original file line number Diff line number Diff line change
Expand Up @@ -198,7 +198,7 @@ def generate_geometry_tables(self):
ufl_geometry = {
ufl.geometry.FacetEdgeVectors: "facet_edge_vertices",
ufl.geometry.CellFacetJacobian: "cell_facet_jacobian",
ufl.geometry.CellEdgeJacobian: "cell_edge_jacobian",
ufl.geometry.CellRidgeJacobian: "cell_ridge_jacobian",
ufl.geometry.ReferenceCellVolume: "reference_cell_volume",
ufl.geometry.ReferenceFacetVolume: "reference_facet_volume",
ufl.geometry.ReferenceCellEdgeVectors: "reference_cell_edge_vectors",
Expand Down
7 changes: 4 additions & 3 deletions ffcx/codegeneration/symbols.py
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ def entity(self, entity_type, restriction):
return self.entity_local_index[0]
elif entity_type == "vertex":
return self.entity_local_index[0]
elif entity_type == "edge":
elif entity_type == "ridge":
return self.entity_local_index[0]
else:
logging.exception(f"Unknown entity_type {entity_type}")
Expand Down Expand Up @@ -138,8 +138,9 @@ def x_component(self, mt):
def J_component(self, mt):
"""Jacobian component."""
# FIXME: Add domain number!
return L.Symbol(format_mt_name(
f"J{mt.expr.ufl_domain().ufl_id()}", mt), dtype=L.DataType.REAL)
return L.Symbol(
format_mt_name(f"J{mt.expr.ufl_domain().ufl_id()}", mt), dtype=L.DataType.REAL
)

def domain_dof_access(self, dof, component, gdim, num_scalar_dofs, restriction):
"""Domain DOF access."""
Expand Down
Loading

0 comments on commit 7a9c0cf

Please sign in to comment.