Implementing Conforming CR #437 #708
Replies: 9 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) |
Beta Was this translation helpful? Give feedback.
-
|
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. |
Beta Was this translation helpful? Give feedback.
-
|
I can also do it if you want, it will be in July after I finish some other work. |
Beta Was this translation helpful? Give feedback.
-
|
Ok. I will test everything locally and then open a PR. We can take it from there. Thanks |
Beta Was this translation helpful? Give feedback.
-
|
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 += 1Partition 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 |
Beta Was this translation helpful? Give feedback.
-
|
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... |
Beta Was this translation helpful? Give feedback.
-
|
I am able to reproduce example 36 now.
The mistake was in the dof ordering above. Changed it to match that in |
Beta Was this translation helpful? Give feedback.
-
|
Fantastic! I'll review the PR as soon as it's ready. |
Beta Was this translation helpful? Give feedback.
-
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 haveElementTetDGandElementTriDGalready 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)
Beta Was this translation helpful? Give feedback.
All reactions