Skip to content

Apparently working BDMCF/E elements#18

Merged
dham merged 28 commits into
masterfrom
bdmc
Sep 20, 2023
Merged

Apparently working BDMCF/E elements#18
dham merged 28 commits into
masterfrom
bdmc

Conversation

@dham

@dham dham commented Aug 8, 2019

Copy link
Copy Markdown
Member

PR for review only. Don't merge as we need to PR into upstream.

This is linked to the Firedrake PR: firedrakeproject/firedrake#1488

@dham dham mentioned this pull request Aug 8, 2019
wence-
wence- previously requested changes May 1, 2020
Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
for j, f in enumerate(polynomials[:, 0]):
T[j, 0, i] = f.evalf(subs=subs)
for j, f in enumerate(polynomials[:, 1]):
T[j, 1, i] = f.evalf(subs=subs)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should use sympy.lambdify so it's not horribly slow like the existing serendipity stuff.

Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 1)))
cur += 2*triangular_number(degree - 1)

formdegree = 1

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only used in the call to the superclass constructor.

Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
return self._degree

def get_nodal_basis(self):
raise NotImplementedError("get_nodal_basis not implemented for bdmce")

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

BDMCE/F ?

Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
raise NotImplementedError("get_dual_set is not implemented for bdmce")

def get_coeffs(self):
raise NotImplementedError("get_coeffs not implemented for bdmce")

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Even though the construction is done on the prime basis directly (for tabulation), it would be good to implement a dual basis given that these elements are interpolatory.

Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
return self.entity_closure_ids

def value_shape(self):
return (2,)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be ref_el.get_dimension() dependent or something, right?

Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
return int(len(self.basis[(0, 0)])/2)


def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid):

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

docstrings. Maybe make a reference to the paper that introduces these elements?

Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
flat_el = flatten_reference_cube(ref_el)
dim = flat_el.get_spatial_dimension()
if dim != 2:
raise Exception("BDMcf_k elements only valid for dimension 2")

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems strange to call something a cube element when it's only valid for squares.

Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
flat_el = flatten_reference_cube(ref_el)
dim = flat_el.get_spatial_dimension()
if dim != 2:
raise Exception("BDMcf_k elements only valid for dimension 2")

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And here.

@wence-

wence- commented Jul 24, 2020

Copy link
Copy Markdown

Can we add fiat level tests for these elements? It appears they don't produce the right results at lowest order.

@wence-

wence- commented Jul 24, 2020

Copy link
Copy Markdown

I am confused by the association of basis functions to topological entities (at least for BDMce) at degree 1.

The FIAT quad is numbered:

     3    
  1-----3
  |     |
0 |     | 1
  |     |
  0-----2
     2
import FIAT
cell = FIAT.reference_element.UFCQuadrilateral()
element = FIAT.brezzi_douglas_marini_cube.BrezziDouglasMariniCubeEdge(cell, 1)
element.entity_dofs()
=> {0: {0: [], 1: [], 2: [], 3: []},
 1: {0: [0, 1], 1: [2, 3], 2: [4, 5], 3: [6, 7]},
 2: {0: []}}

So let's look at the first two basis functions which lie on the edge x = 0.

print(element.basis[0, 0][:2])
=> [[0, 1.0*x - 1.0], [0.5*y*(1.0*y - 1.0), (1.0 - 1.0*x)*(1.0 - 2*y)]]

I don't know exactly what these should be, but I was expecting things that have a x component that varies linearly in y, and are constant in x. Is this wrong?

In any case, if I look at the next two basis functions

print(element.basis[0,0][2:4])                                         
=> [[0, -1.0*x], [0.5*y*(1.0 - 1.0*y), 1.0*x*(1.0 - 2*y)]]

It seems these pair up "nicely?" with the corresponding basis functions 0 and 1. What am I missing?

Given this, it would be good if the implementation pointed to a reference for the basis functions. If they're the standard ones described in the Brezzi Boffi and Fortin book, it seems odd that they're defined in this primal way, rather than by using a dual basis (which is explicitly described there).

Comment thread FIAT/brezzi_douglas_marini_cube.py Outdated
[(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))] +
[(-leg(j, x_mid)*dy[1], 0) for j in range(deg)] +
[(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))])

@tommbendall tommbendall Jan 25, 2021

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My understanding of the issues @pefarrell was reporting with the lowest-order spaces is that these basis functions aren't quite right.

When taking the 2D curl of each basis function (or the divergence for BDMCF) it should return a polynomial of deg-1.

Inspection of these basis functions shows that this is correct for the functions such as
[(0, -leg(j, y_mid)*dx[0]) for j in range(deg)]
but not quite for those like
[(-leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[0])].

I did a bit of maths and I think that one component of these needs modifying by a coefficient. Running the bdmc_riesz test in the bdmc branch of Firedrake, I got a convergence rate of deg+1 when using the following functions (which ensure that the curl/div of the basis functions is of the correct polynomial order):

from sympy import binomial
coeff = binomial(2*deg,deg) / ((deg+1)*binomial(2*(deg-1), deg-1))
EL = tuple([(0, -leg(j, y_mid)*dx[0]) for j in range(deg)] +
           [(coeff*-leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[0])] +
           [(0, -leg(j, y_mid)*dx[1]) for j in range(deg)] +
           [(coeff*leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[1])] +
           [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] +
           [(-leg(deg, x_mid)*dy[0], coeff*-leg(deg-1, x_mid)*dx[0]*dx[1])] +
           [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)] +
           [(-leg(deg, x_mid)*dy[1], coeff*leg(deg-1, x_mid)*dx[0]*dx[1])])

@tommbendall

tommbendall commented Jan 28, 2021

Copy link
Copy Markdown

I think I've addressed most of @wence- 's comments -- I have not changed the name cube (even though these are elements on quads), and I have not yet implemented a dual basis.

@jmv2009

jmv2009 commented Jun 6, 2021

Copy link
Copy Markdown

FYI: Just is in #22, this will not work in Fenics without dual basis.

@wence-

wence- commented Jun 29, 2022

Copy link
Copy Markdown

This was superseded by #22 and can be closed, yes? Ping @rckirby

@rckirby

rckirby commented Jun 29, 2022

Copy link
Copy Markdown

My understanding is that BDMCE/F implement the S elements, while trimmed serendipity implement the S- elements so that they are a bit different. Maybe @dham recalls what was going on there?

@dham
dham merged commit e440a06 into master Sep 20, 2023
@dham
dham deleted the bdmc branch September 20, 2023 15:15
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

7 participants