forked from FEniCS/fiat
-
Notifications
You must be signed in to change notification settings - Fork 7
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
67 changed files
with
2,027 additions
and
930 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,84 @@ | ||
| # Copyright (C) 2024 Pablo D. Brubeck | ||
| # | ||
| # This file is part of FIAT (https://www.fenicsproject.org) | ||
| # | ||
| # SPDX-License-Identifier: LGPL-3.0-or-later | ||
| # | ||
| # Written by Pablo D. Brubeck ([email protected]), 2024 | ||
|
|
||
| from FIAT import finite_element, dual_set, polynomial_set | ||
| from FIAT.functional import ComponentPointEvaluation, PointDivergence | ||
| from FIAT.quadrature_schemes import create_quadrature | ||
| from FIAT.macro import CkPolynomialSet, AlfeldSplit | ||
|
|
||
| import numpy | ||
|
|
||
|
|
||
| def AlfeldSorokinaSpace(ref_el, degree): | ||
| """Return a vector-valued C0 PolynomialSet on an Alfeld split with C0 | ||
| divergence. This works on any simplex and for all polynomial degrees.""" | ||
| ref_complex = AlfeldSplit(ref_el) | ||
| sd = ref_complex.get_spatial_dimension() | ||
| C0 = CkPolynomialSet(ref_complex, degree, order=0, shape=(sd,), variant="bubble") | ||
| expansion_set = C0.get_expansion_set() | ||
| num_members = C0.get_num_members() | ||
| coeffs = C0.get_coeffs() | ||
|
|
||
| facet_el = ref_complex.construct_subelement(sd-1) | ||
| phi = polynomial_set.ONPolynomialSet(facet_el, 0 if sd == 1 else degree-1) | ||
| Q = create_quadrature(facet_el, 2 * phi.degree) | ||
| qpts, qwts = Q.get_points(), Q.get_weights() | ||
| phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] | ||
| weights = numpy.multiply(phi_at_qpts, qwts) | ||
|
|
||
| rows = [] | ||
| for facet in ref_complex.get_interior_facets(sd-1): | ||
| n = ref_complex.compute_normal(facet) | ||
| jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=1) | ||
| div_jump = n[:, None, None] * jumps[1][None, ...] | ||
| r = numpy.tensordot(div_jump, weights, axes=(-1, -1)) | ||
| rows.append(r.reshape(num_members, -1).T) | ||
|
|
||
| if len(rows) > 0: | ||
| dual_mat = numpy.vstack(rows) | ||
| nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) | ||
| coeffs = numpy.tensordot(nsp, coeffs, axes=(-1, 0)) | ||
| return polynomial_set.PolynomialSet(ref_complex, degree, degree, expansion_set, coeffs) | ||
|
|
||
|
|
||
| class AlfeldSorokinaDualSet(dual_set.DualSet): | ||
| def __init__(self, ref_el, degree): | ||
| if degree != 2: | ||
| raise NotImplementedError(f"{type(self).__name__} only defined for degree = 2") | ||
|
|
||
| top = ref_el.get_topology() | ||
| sd = ref_el.get_spatial_dimension() | ||
| entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} | ||
|
|
||
| nodes = [] | ||
| for dim in sorted(top): | ||
| for entity in sorted(top[dim]): | ||
| cur = len(nodes) | ||
|
|
||
| dpts = ref_el.make_points(dim, entity, degree-1) | ||
| nodes.extend(PointDivergence(ref_el, pt) for pt in dpts) | ||
|
|
||
| pts = ref_el.make_points(dim, entity, degree) | ||
| nodes.extend(ComponentPointEvaluation(ref_el, k, (sd,), pt) | ||
| for pt in pts for k in range(sd)) | ||
| entity_ids[dim][entity].extend(range(cur, len(nodes))) | ||
|
|
||
| super().__init__(nodes, ref_el, entity_ids) | ||
|
|
||
|
|
||
| class AlfeldSorokina(finite_element.CiarletElement): | ||
| """The Alfeld-Sorokina C0 quadratic macroelement with C0 divergence. | ||
| This element belongs to a Stokes complex, and is paired with CG1(Alfeld). | ||
| """ | ||
| def __init__(self, ref_el, degree=2): | ||
| dual = AlfeldSorokinaDualSet(ref_el, degree) | ||
| poly_set = AlfeldSorokinaSpace(ref_el, degree) | ||
| formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form | ||
| super().__init__(poly_set, dual, degree, formdegree, | ||
| mapping="contravariant piola") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,71 @@ | ||
| # This file is part of FIAT (https://www.fenicsproject.org) | ||
| # | ||
| # SPDX-License-Identifier: LGPL-3.0-or-later | ||
| # | ||
| # Written by Pablo D. Brubeck ([email protected]), 2024 | ||
|
|
||
| from FIAT import finite_element, polynomial_set | ||
| from FIAT.bernardi_raugel import BernardiRaugelDualSet | ||
| from FIAT.quadrature_schemes import create_quadrature | ||
| from FIAT.reference_element import TRIANGLE | ||
| from FIAT.macro import CkPolynomialSet | ||
| from FIAT.hct import HsiehCloughTocher | ||
|
|
||
| import numpy | ||
|
|
||
|
|
||
| def ArnoldQinSpace(ref_el, degree, reduced=False): | ||
| """Return a basis for the Arnold-Qin space | ||
| curl(HCT-red) + P0 x if reduced = True, and | ||
| curl(HCT) + P0 x if reduced = False.""" | ||
| if ref_el.get_shape() != TRIANGLE: | ||
| raise ValueError("Arnold-Qin only defined on triangles") | ||
| if degree != 2: | ||
| raise ValueError("Arnold-Qin only defined for degree = 2") | ||
| sd = ref_el.get_spatial_dimension() | ||
| HCT = HsiehCloughTocher(ref_el, degree+1, reduced=True) | ||
| ref_complex = HCT.get_reference_complex() | ||
| Q = create_quadrature(ref_complex, 2 * degree) | ||
| qpts, qwts = Q.get_points(), Q.get_weights() | ||
|
|
||
| x = qpts.T | ||
| bary = numpy.asarray(ref_el.make_points(sd, 0, sd+1)) | ||
| P0x_at_qpts = x[None, :, :] - bary[:, :, None] | ||
|
|
||
| tab = HCT.tabulate(1, qpts) | ||
| curl_at_qpts = numpy.stack([tab[(0, 1)], -tab[(1, 0)]], axis=1) | ||
| if reduced: | ||
| curl_at_qpts = curl_at_qpts[:9] | ||
|
|
||
| C0 = CkPolynomialSet(ref_complex, degree, order=0, scale=1, variant="bubble") | ||
| C0_at_qpts = C0.tabulate(qpts)[(0,) * sd] | ||
| duals = numpy.multiply(C0_at_qpts, qwts) | ||
| M = numpy.dot(duals, C0_at_qpts.T) | ||
| duals = numpy.linalg.solve(M, duals) | ||
|
|
||
| # Remove the constant nullspace | ||
| ids = [0, 3, 6] | ||
| A = numpy.asarray([[1, 1, 1], [1, -1, 0], [0, -1, 1]]) | ||
| phis = curl_at_qpts | ||
| phis[ids] = numpy.tensordot(A, phis[ids], axes=(-1, 0)) | ||
| # Replace the constant nullspace with P_0 x | ||
| phis[0] = P0x_at_qpts | ||
| coeffs = numpy.tensordot(phis, duals, axes=(-1, -1)) | ||
| return polynomial_set.PolynomialSet(ref_complex, degree, degree, | ||
| C0.get_expansion_set(), coeffs) | ||
|
|
||
|
|
||
| class ArnoldQin(finite_element.CiarletElement): | ||
| """The Arnold-Qin C0(Alfeld) quadratic macroelement with divergence in P0. | ||
| This element belongs to a Stokes complex, and is paired with unsplit DG0.""" | ||
| def __init__(self, ref_el, degree=2, reduced=False): | ||
| poly_set = ArnoldQinSpace(ref_el, degree) | ||
| if reduced: | ||
| order = 1 | ||
| mapping = "contravariant piola" | ||
| else: | ||
| order = degree | ||
| mapping = "affine" | ||
| dual = BernardiRaugelDualSet(ref_el, order, degree=degree) | ||
| formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form | ||
| super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) |
Oops, something went wrong.