firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
496 stars 157 forks source link

Making a broken element on an extruded sphere mesh #1762

Closed colinjcotter closed 4 years ago

colinjcotter commented 4 years ago

The following code fails an exception when trying to get the quadrature degree when making a broken element. What went wrong?

from firedrake import *
import ufl
Radius = 100
layers = 1
layer_height = 1
refinement = 4
meshdegree = 1
hdegree = 0
vdegree = 0
m = IcosahedralSphereMesh(radius=Radius, refinement_level=refinement, degree=meshdegree)
mesh = ExtrudedMesh(m, layers, layer_height=layer_height, extrusion_type='radial')
BDM = FiniteElement("BDM", "triangle", hdegree+1)
DP_ = FiniteElement("DP", "interval", vdegree+0)
HDiv1 = HDiv(TensorProductElement(BDM, DP_))
DP = FiniteElement("DP", "triangle", hdegree+0)
P_ = FiniteElement("P", "interval", vdegree+1)
HDiv2 = HDiv(TensorProductElement(DP, P_))
W2element = EnrichedElement(HDiv1 + HDiv2)
W2 = FunctionSpace(mesh, W2element)
U0 = Function(W2) #velocity
brokenU = ufl.BrokenElement(W2element)
V_d = FunctionSpace(mesh, brokenU)
wence- commented 4 years ago

PEBKAC, I think, write W2element = EnrichedElement(HDiv1, HDiv2) or W2element = HDiv1 + HDiv2 (which both make an enriched element, not W2element = EnrichedElement(HDiv1 + HDiv2)

wence- commented 4 years ago

This patch in finat fixes this error by unwinding enrichedelement with only one argument at construction time:

diff --git a/finat/enriched.py b/finat/enriched.py
index c494029..da0e7d1 100644
--- a/finat/enriched.py
+++ b/finat/enriched.py
@@ -13,9 +13,13 @@ class EnrichedElement(FiniteElementBase):
     """A finite element whose basis functions are the union of the
     basis functions of several other finite elements."""

-    def __init__(self, elements):
-        super(EnrichedElement, self).__init__()
-        self.elements = tuple(elements)
+    def __new__(cls, elements):
+        if len(elements) == 1:
+            return elements[0]
+        else:
+            self = super().__new__(cls)
+            self.elements = tuple(elements)
+            return self

     @cached_property
     def cell(self):
colinjcotter commented 4 years ago

Ah, thanks! I'll tell my project student.


From: Lawrence Mitchell notifications@github.com Sent: 01 July 2020 17:46 To: firedrakeproject/firedrake firedrake@noreply.github.com Cc: Cotter, Colin J colin.cotter@imperial.ac.uk; Author author@noreply.github.com Subject: Re: [firedrakeproject/firedrake] Making a broken element on an extruded sphere mesh (#1762)

PEBKAC, I think, write W2element = EnrichedElement(HDiv1, HDiv2) or W2element = HDiv1 + HDiv2 (which both make an enriched element, not W2element = EnrichedElement(HDiv1 + HDiv2)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/firedrakeproject/firedrake/issues/1762#issuecomment-652530117, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABOSV4WPGQA4JIDJHHG27H3RZNR4XANCNFSM4OLODYMQ.