Replies: 11 comments
-
|
I could not write testcode and copy the codes for 15-parameter Triangle Plate Element as follow. and I want to know whether it is correct. class ElementTri15ParameterPlate(ElementGlobal): |
Beta Was this translation helpful? Give feedback.
-
|
class ElementTriMZ(ElementGlobal): The codes above is for modified Zienkiewicz element, but it doesn't work. |
Beta Was this translation helpful? Give feedback.
-
|
class ElementTriMorleyZienkiewicz(ElementGlobal): The codes above is for Morley-Zienkiewicz element, but it also doesn't work. |
Beta Was this translation helpful? Give feedback.
-
|
I'll take a look if I can figure out what is the issue. |
Beta Was this translation helpful? Give feedback.
-
|
I changed the DOF order in 15 parameter plate element and now it seems to work better: from skfem import *
from skfem.element.element_global import ElementGlobal
from skfem.refdom import RefTri
from skfem.helpers import *
import numpy as np
class ElementTri15ParameterPlate(ElementGlobal):
"""15-parameter Triangle Plate Element"""
nodal_dofs = 3
facet_dofs = 2
maxdeg = 4
dofnames = ['u', 'u_x', 'u_y', 'u', 'u_n']
doflocs = np.array([[0., 0.],
[0., 0.],
[0., 0.],
[1., 0.],
[1., 0.],
[1., 0.],
[0., 1.],
[0., 1.],
[0., 1.],
[.5, 0.],
[.5, .5],
[0., .5],
[.5, 0.],
[.5, .5],
[0., .5]])
refdom = RefTri
def gdof(self, F, w, i):
if i == 0:
return F[()](*w['v'][0])
elif i == 1:
return F[(0,)](*w['v'][0])
elif i == 2:
return F[(1,)](*w['v'][0])
elif i == 3:
return F[()](*w['v'][1])
elif i == 4:
return F[(0,)](*w['v'][1])
elif i == 5:
return F[(1,)](*w['v'][1])
elif i == 6:
return F[()](*w['v'][2])
elif i == 7:
return F[(0,)](*w['v'][2])
elif i == 8:
return F[(1,)](*w['v'][2])
elif i == 9:
return F[()](*w['e'][0])
elif i == 10:
return (F[(0,)](*w['e'][0]) * w['n'][0, 0]
+ F[(1,)](*w['e'][0]) * w['n'][0, 1])
elif i == 11:
return F[()](*w['e'][1])
elif i == 12:
return (F[(0,)](*w['e'][1]) * w['n'][1, 0]
+ F[(1,)](*w['e'][1]) * w['n'][1, 1])
elif i == 13:
return F[()](*w['e'][2])
elif i == 14:
return (F[(0,)](*w['e'][2]) * w['n'][2, 0]
+ F[(1,)](*w['e'][2]) * w['n'][2, 1])
self._index_error()
m = MeshTri().refined(3)
ib = Basis(m, ElementTri15ParameterPlate())
t = 1.
E = 1.
nu = 0.3
D = E * t ** 3 / (12. * (1. - nu ** 2))
@BilinearForm
def bilinf(u, v, w):
def C(T):
trT = T[0, 0] + T[1, 1]
return E / (1. + nu) * \
np.array([[T[0, 0] + nu / (1. - nu) * trT, T[0, 1]],
[T[1, 0], T[1, 1] + nu / (1. - nu) * trT]])
return t ** 3 / 12.0 * ddot(C(dd(u)), dd(v))
def load(x):
return np.sin(np.pi * x[0]) * np.sin(np.pi * x[1])
@LinearForm
def linf(v, w):
return load(w.x) * v
K = asm(bilinf, ib)
f = asm(linf, ib)
x = solve(*condense(K, f, D=ib.get_dofs().all('u')))
from skfem.visuals.matplotlib import *
plot(ib, x, Nrefs=3)
show()
|
Beta Was this translation helpful? Give feedback.
-
|
I think the other methods use a different type of polynomial basis. E.g., Bell triangle does not have full polynomial space and I have not looked into how ElementGlobal should be modified. |
Beta Was this translation helpful? Give feedback.
-
|
Should we add 15 parameter nonconforming triangle to scikit-fem? |
Beta Was this translation helpful? Give feedback.
-
|
Sure, we have more choices when solving four order problems |
Beta Was this translation helpful? Give feedback.
-
|
For the modified Zienkiewicz element, a H2-nonconforming elment, can't be worked by using ElementGlobal class. And now I do another try. I write it in explicit basis and modify ElementH1 class to get a new ElementH2 class, with adding a hessian operator. import numpy as np from ...refdom import RefTri class ElementTriMZ(ElementH2): |
Beta Was this translation helpful? Give feedback.
-
|
I am not sure whether this code is correct. And I compared Hermite element using the same way with ElementTriHermite in scikit-fem. I copy the code of ElementTriHermite with explicit basis here: import numpy as np from ..element_h2 import ElementH2 class ElementTriHermite0(ElementH2): |
Beta Was this translation helpful? Give feedback.
-
|
And they are different to solve the same problem, could you help me look at it and check it? Thank you in advance. |
Beta Was this translation helpful? Give feedback.
-
I find that ElementGlobal is a powerful class. I have ever programmed Argyris element in MATLAB, it is a tedious project. By contrast, ElementTriArgyris class in scikit-fem is concise and readable, it is superb job. For fourth order problem, there are Morley element, Hermite element, and Argyris element over trianglular elements using ElementGlobal. So I take a try for more elements, from the following book,
"Z. Shi, M. Wang, Finite Element Methods, Scientific Publishers, Beijing, 2013."
I find these H2-nonconforming elements for four order problem. they are Zienkiewicz element, Morley-Zienkiewicz element, Modified Zienkiewicz element, 12-parameter Triangle Plate Element and 15-parameter Triangle Plate Element, but only last one works.
I'm wondering why the others cannot be implemented by using ElementGlobal. I guess the reason is that the number of DOFs of Morley elment, Hermite element, 15-parameter Triangle Plate element and Argyris element is 6,10,15,21, corresponding to {1,x,y,x^2,xy,y^2},{1,x,y,x^2,xy,y^2,x^3,x^2y,xy^2,y^3},{1,x,y,x^2,xy,y^2,x^3,x^2y,xy^2,y^3,x^4,x^3y,x^2y^2,xy^3,y^4} and {1,x,y,x^2,xy,y^2,x^3,x^2y,xy^2,y^3,x^4,x^3y,x^2y^2,xy^3,y^4,x^5,x^4y,x^3y^2,x^2y^3,xy^4,y^5}.
My question is how to implement the rest elements.
Beta Was this translation helpful? Give feedback.
All reactions