diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py index 2c49ddd8967f63013cdd3f4c4ec5cdfe9918196c..66d3179d91e34e670ef8adfd6db90b3fe19d4459 100644 --- a/generate_all_hyteg_forms.py +++ b/generate_all_hyteg_forms.py @@ -33,6 +33,7 @@ from hog.element_geometry import ( ) from hog.function_space import ( LagrangianFunctionSpace, + DGFunctionSpace, N1E1Space, P2PlusBubbleSpace, TrialSpace, @@ -132,14 +133,34 @@ class FormInfo: def space_desc(self) -> str: """A compact representation of the function spaces.""" + descr_string = "" if self.trial_family == "N1E1": - return "n1e1" + descr_string = "n1e1" elif self.trial_family == "P2 enhanced with Bubble": - return "p2_plus_bubble" - elif self.trial_degree == self.test_degree: - return f"p{self.trial_degree}" - else: - return f"p{self.trial_degree}_to_p{self.test_degree}" + if self.test_family == "P2 enhanced with Bubble": + descr_string = "p2_plus_bubble" + elif self.test_family == "DG": + descr_string = f"p2_plus_bubble_to_dg{self.test_degree}" + elif self.trial_family == "Lagrange" and self.test_family == "Lagrange": + if self.trial_degree == self.test_degree: + descr_string = f"p{self.trial_degree}" + else: + descr_string = f"p{self.trial_degree}_to_p{self.test_degree}" + elif self.trial_family == "DG": + if self.test_family == "DG": + if self.trial_degree == self.test_degree: + descr_string = f"dg{self.trial_degree}" + else: + descr_string = f"dg{self.trial_degree}_to_dg{self.test_degree}" + elif self.test_family == "Lagrange": + descr_string = f"dg{self.trial_degree}_to_p{self.test_degree}" + elif self.test_family == "P2 enhanced with Bubble": + descr_string = f"dg{self.trial_degree}_to_p2_plus_bubble" + else: + raise HOGException( + f"Do not know how to name combination of DGFunctionSpace with {self.test_family}." + ) + return descr_string def blending_desc(self) -> str: """The type of transformation from the reference element. Either 'blending' or 'affine'.""" @@ -634,6 +655,48 @@ for trial_deg, test_deg, transpose in [ ) ) +for trial_deg, test_deg, transpose in [ + (1, 2, True), + (2, 1, False), +]: + for blending in [IdentityMap(), ExternalMap()]: + if not transpose: + form_infos.append( + FormInfo( + f"div", + trial_degree=trial_deg, + test_degree=test_deg, + trial_family="P2 enhanced with Bubble", + test_family="DG", + quad_schemes={ + 2: 3 + test_deg - 1, + 3: 4 + test_deg - 1, + }, + row_dim=1, + col_dim=3, + is_implemented=is_implemented_for_vector_to_scalar, + blending=blending, + ) + ) + else: + form_infos.append( + FormInfo( + f"divt", + trial_degree=trial_deg, + test_degree=test_deg, + trial_family="DG", + test_family="P2 enhanced with Bubble", + quad_schemes={ + 2: trial_deg + 3 - 1, + 3: trial_deg + 4 - 1, + }, + row_dim=3, + col_dim=1, + is_implemented=is_implemented_for_scalar_to_vector, + blending=blending, + ) + ) + for blending in [IdentityMap(), ExternalMap()]: form_infos.append( FormInfo( diff --git a/hog/function_space.py b/hog/function_space.py index 930e200b0fd6430575bd12b49099fd7e09ae611c..2bd124b1fc5d26764f888beee15133a323306c74 100644 --- a/hog/function_space.py +++ b/hog/function_space.py @@ -699,3 +699,166 @@ class P2PlusBubbleSpace(FunctionSpace): def __repr__(self): return str(self) + + +class DGFunctionSpace(FunctionSpace): + """ + Space of discontinuous piecewise polynomial functions (Discontinuous-Galerkin) + """ + + def __init__(self, degree: int, symbolizer: Symbolizer): + """Creates a FunctionSpace of the Discontinuous Galerkin family. + + :param degree: the order of the shape functions + :param symbolizer: a Symbolizer instance + """ + if degree not in [1, 2]: + raise HOGException("Only degrees 1 and 2 are supported.") + + self._degree = degree + self._symbolizer = symbolizer + + @property + def is_vectorial(self) -> bool: + return False + + @property + def is_continuous(self) -> bool: + return False + + @property + def degree(self) -> int: + return self._degree + + @property + def symbolizer(self) -> Symbolizer: + return self._symbolizer + + @property + def family(self) -> str: + return "DG" + + def shape( + self, + geometry: ElementGeometry, + domain: str = "reference", + dof_map: Optional[List[int]] = None, + ) -> List[sp.MatrixBase]: + """Returns a list containing the shape functions on the element. + + :param dof_map: this list can be used to specify (remap) the DoF ordering of the element + """ + if domain in ["ref", "reference"]: + symbols = self.symbolizer.ref_coords_as_list(geometry.dimensions) + + basis_functions = [] + + if ( + isinstance(geometry, TriangleElement) + and self.family in ["DG"] + and self._degree == 1 + ): + basis_functions = [ + 1 - symbols[0] - symbols[1], + symbols[0], + symbols[1], + ] + + elif ( + isinstance(geometry, TriangleElement) + and self.family in ["DG"] + and self._degree == 2 + ): + x = symbols[0] + y = symbols[1] + basis_functions = [ + 2 * x**2 + 4 * x * y - 3 * x + 2 * y**2 - 3 * y + 1, + 2 * x**2 - x, + 2 * y**2 - y, + 4 * x * y, + -4 * x * y - 4 * y**2 + 4 * y, + -4 * x**2 - 4 * x * y + 4 * x, + ] + + elif ( + isinstance(geometry, TetrahedronElement) + and self.family in ["DG"] + and self._degree == 1 + ): + basis_functions = [ + 1 - symbols[0] - symbols[1] - symbols[2], + symbols[0], + symbols[1], + symbols[2], + ] + + elif ( + isinstance(geometry, TetrahedronElement) + and self.family in ["DG"] + and self._degree == 2 + ): + xi_1 = symbols[0] + xi_2 = symbols[1] + xi_3 = symbols[2] + basis_functions = [ + ( + 2.0 * xi_1 * xi_1 + + 4.0 * xi_1 * xi_2 + + 4.0 * xi_1 * xi_3 + - 3.0 * xi_1 + + 2.0 * xi_2 * xi_2 + + 4.0 * xi_2 * xi_3 + - 3.0 * xi_2 + + 2.0 * xi_3 * xi_3 + - 3.0 * xi_3 + + 1.0 + ), + (2.0 * xi_1 * xi_1 - 1.0 * xi_1), + (2.0 * xi_2 * xi_2 - 1.0 * xi_2), + (2.0 * xi_3 * xi_3 - 1.0 * xi_3), + (4.0 * xi_2 * xi_3), + (4.0 * xi_1 * xi_3), + (4.0 * xi_1 * xi_2), + ( + -4.0 * xi_1 * xi_3 + - 4.0 * xi_2 * xi_3 + - 4.0 * xi_3 * xi_3 + + 4.0 * xi_3 + ), + ( + -4.0 * xi_1 * xi_2 + - 4.0 * xi_2 * xi_2 + - 4.0 * xi_2 * xi_3 + + 4.0 * xi_2 + ), + ( + -4.0 * xi_1 * xi_1 + - 4.0 * xi_1 * xi_2 + - 4.0 * xi_1 * xi_3 + + 4.0 * xi_1 + ), + ] + + else: + raise HOGException( + "Basis functions not implemented for the specified element type and geometry." + ) + + if dof_map: + raise HOGException("DoF reordering not implemented.") + + basis_functions = [sp.Matrix([b]) for b in basis_functions] + return basis_functions + + raise HOGException(f"Shape function not available for domain type {domain}") + + def __eq__(self, other: Any) -> bool: + if not isinstance(other, LagrangianFunctionSpace): + return False + return self.family == other.family and self._degree == other._degree + + def __str__(self): + return f"{self.family}, degree: {self._degree}" + + def __repr__(self): + return str(self)