-
Notifications
You must be signed in to change notification settings - Fork 80
This issue was moved to a discussion.
You can continue the conversation there. Go to discussion →
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implementing Conforming CR #437 #651
Comments
The basis functions are copied from (124) in this paper which I have verified to be correct by implementing it in a commercial code (as Abaqus user subroutine to be specific) |
Here are some tests that are run against the elements (if applicable):
We also test the convergence rates, e.g.: https://github.com/kinnala/scikit-fem/blob/master/tests/test_convergence.py#L213 A new file is created: If you want you can open a pull request with these changes. |
I can also do it if you want, it will be in July after I finish some other work. |
Ok. I will test everything locally and then open a PR. We can take it from there. Thanks |
Ok some progress... Kronecker DeltaIn [1]: from skfem.element.element_tet import ElementTetCCR as tetCR
In [2]: elemtet = tetCR()
In [3]: import numpy as np
...: doflocs = elemtet.doflocs
...: num_dofs = doflocs.shape[0]
...: basis_at_doflocs = np.array([elemtet.lbasis(doflocs.T, i)[0] for i in range(num_dofs)], float)
In [4]: np.testing.assert_allclose(basis_at_doflocs, np.eye(num_dofs), atol=1.e-14) Derivative valuesIn [5]: eps=1.e-6
In [6]: for base in [0., .3, .6, .9]:
...: if elemtet.dim == 1:
...: y = np.array([[base, base + eps]])
...: elif elemtet.dim == 2:
...: y = np.array([[base, base + eps, base, base],
...: [base, base, base, base + eps]])
...: elif elemtet.dim == 3:
...: y = np.array([[base, base + eps, base, base, base, base],
...: [base, base, base, base + eps, base, base],
...: [base, base, base, base, base, base + eps]])
...: i = 0
...: while True:
...: try:
...: out = elemtet.lbasis(y, i)
...: except ValueError:
...: break
...: diff = (out[0][1] - out[0][0]) / eps
...: errmsg = 'x-derivative for {}th bfun failed for {}'
...: np.testing.assert_almost_equal(diff, out[1][0][0], decimal=3,
...: err_msg=errmsg.format(i, elemtet))
...: if elemtet.dim > 1:
...: diff = (out[0][3] - out[0][2]) / eps
...: errmsg = 'y-derivative for {}th bfun failed for {}'
...: np.testing.assert_almost_equal(diff, out[1][1][3], decimal=3,
...: err_msg=errmsg.format(i, elemtet))
...: if elemtet.dim == 3:
...: diff = (out[0][5] - out[0][4]) / eps
...: errmsg = 'z-derivative for {}th bfun failed for {}'
...: np.testing.assert_almost_equal(diff, out[1][2][4], decimal=3,
...: err_msg=errmsg.format(i, elemtet))
...: i += 1 Partition of UnityIn [7]: if elemtet.dim == 1:
...: y = np.array([[.15]])
...: elif elemtet.dim == 2:
...: y = np.array([[.15],
...: [.15]])
...: elif elemtet.dim == 3:
...: y = np.array([[.15],
...: [.15],
...: [.15]])
...: out = 0.
...: for i in range(elemtet.doflocs.shape[0]):
...: out += elemtet.lbasis(y, i)[0][0]
...: np.testing.assert_almost_equal(out, 1, err_msg='failed for {}'.format(elemtet)) and similarly for the |
There are still some persistent issues (for e.g., Newton not converging for larger values of the applied stretches, whereas it does in the case of the Taylor-Hood approximation). I used a higher degree quadrature rule, namely quad_points = np.array([[0.25 , 0. , 0.33333333, 0.33333333, 0.33333333,
0.72727273, 0.09090909, 0.09090909, 0.09090909, 0.43344985,
0.06655015, 0.06655015, 0.06655015, 0.43344985, 0.43344985],
[0.25 , 0.33333333, 0.33333333, 0.33333333, 0. ,
0.09090909, 0.09090909, 0.09090909, 0.72727273, 0.06655015,
0.43344985, 0.06655015, 0.43344985, 0.06655015, 0.43344985],
[0.25 , 0.33333333, 0.33333333, 0. , 0.33333333,
0.09090909, 0.09090909, 0.72727273, 0.09090909, 0.06655015,
0.06655015, 0.43344985, 0.43344985, 0.43344985, 0.06655015]
])
quad_weights = np.array([0.03028368, 0.00602679, 0.00602679, 0.00602679, 0.00602679,
0.01164525, 0.01164525, 0.01164525, 0.01164525, 0.01094914,
0.01094914, 0.01094914, 0.01094914, 0.01094914, 0.01094914])
mesh = MeshTet().refined(2)
uelem = ElementVectorH1(ElementTetCCR())
pelem = ElementTetDG(ElementTetP1())
....
....
basis = {
field: InteriorBasis(mesh, e, quadrature=(quad_points, quad_weights))
for field, e in elems.items()
} Will investigate further... |
Fantastic! I'll review the PR as soon as it's ready. |
Merged |
This issue was moved to a discussion.
You can continue the conversation there. Go to discussion →
Hi @kinnala,
I finally got back to implementing conforming CR (see 3.2.3 in here for a concise definition of the spaces for velocity and pressure) using the suggestions in #437 for incompressible stokes (hyperelasticity). Since I plan to use this element in my work, and if you deem it fit then include it in
scikit-fem
. This is what I have for the definitions so far...3D
2D
I am looking to first test this in example 36 (in doing so I realized that example 36 was written in a non-standard fashion -- I plan to clean that up soon as well). Since the pressure space is discontinuous
P1
, and now that you haveElementTetDG
andElementTriDG
already inscikit-fem
, it seems that example 36 should be reproducible by simply changing the relevant spaces toPlease let me know what you think. In the meanwhile, I will check if I have implemented everything correctly (basis functions and their gradients)
The text was updated successfully, but these errors were encountered: