bhaveshshrimali commented 3 years ago

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

class ElementTetCCR(ElementH1):

    nodal_dofs = 1
    facet_dofs = 1
    edge_dofs = 1
    interior_dofs = 1
    maxdeg = 4
    dofnames = ["u", "u", "u", "u"]
    doflocs = np.array(
        [
            [1.0, 0.0, 0.0],
            [0.0, 1.0, 0.0],
            [0.0, 0.0, 0.0],
            [0.0, 0.0, 1.0],
            [0.5, 0.5, 0.0],
            [0.0, 0.5, 0.0],
            [0.5, 0.0, 0.0],
            [0.5, 0.0, 0.5],
            [0.0, 0.5, 0.5],
            [0.0, 0.0, 0.5],
            [1.0 / 3, 1.0 / 3, 0.0],
            [1.0 / 3, 1.0 / 3, 1.0 / 3],
            [0.0, 1.0 / 3, 1.0 / 3],
            [1.0 / 3, 0.0, 1.0 / 3],
            [1.0 / 4, 1.0 / 4, 1.0 / 4],
        ]
    )
    refdom = RefTet

    def lbasis(self, X, i):
        x, y, z = X 

        if i == 0:  # at (1,0,0)
            phi = 3 * x * (y * z + (y + z) * (-x - y - z + 1)) + x * (
                2 * x - 4 * y * z * (-x - y - z + 1) - 1
            )
            dphi = np.array(
                [
                    -3*x*(y + z) + 2*x*(2*y*z + 1) + 2*x + 4*y*z*(x + y + z - 1) + 3*y*z - 3*(y + z)*(x + y + z - 1) - 1,
                    x*(4*z - 3)*(x + 2*y + z - 1),
                    x*(4*y - 3)*(x + y + 2*z - 1),
                ]
            )
        elif i == 1:  # at (1,0,0)
            phi = y*(4*x*z*(x + y + z - 1) + 3*x*z + 2*y - 3*(x + z)*(x + y + z - 1) - 1)
            dphi = np.array(
                [
                    y*(4*z - 3)*(2*x + y + z - 1),
                    4*x*z*(x + y + z - 1) + 3*x*z - 3*y*(x + z) + 2*y*(2*x*z + 1) + 2*y - 3*(x + z)*(x + y + z - 1) - 1,
                    y*(4*x - 3)*(x + y + 2*z - 1),
                ]
            )
        elif i == 2:  # at (0,1,0)
            phi = (x + y + z - 1)*(4*x*y*z - 3*x*y - 3*x*z + 2*x - 3*y*z + 2*y + 2*z - 1)
            dphi = np.array(
                [
                    8*x*y*z - 6*x*y - 6*x*z + 4*x + 4*y**2*z - 3*y**2 + 4*y*z**2 - 13*y*z + 7*y - 3*z**2 + 7*z - 3,
                    4*x**2*z - 3*x**2 + 8*x*y*z - 6*x*y + 4*x*z**2 - 13*x*z + 7*x - 6*y*z + 4*y - 3*z**2 + 7*z - 3,
                    4*x**2*y - 3*x**2 + 4*x*y**2 + 8*x*y*z - 13*x*y - 6*x*z + 7*x - 3*y**2 - 6*y*z + 7*y + 4*z - 3,
                ]
            )
        elif i == 3:  # at (0,0,1)
            phi = z*(4*x*y*(x + y + z - 1) + 3*x*y + 2*z - 3*(x + y)*(x + y + z - 1) - 1)
            dphi = np.array(
                [
                    z*(4*y - 3)*(2*x + y + z - 1),
                    z*(4*x - 3)*(x + 2*y + z - 1),
                    4*x*y*(x + y + z - 1) + 3*x*y - 3*z*(x + y) + 2*z*(2*x*y + 1) + 2*z - 3*(x + y)*(x + y + z - 1) - 1,
                ]
            )
        elif i == 4:  # between (0,1)
            phi = 4*x*y*(3*x + 3*y - 8*z*(x + y + z - 1) - 2)
            dphi = np.array(
                [
                    4*y*(-x*(8*z - 3) + 3*x + 3*y - 8*z*(x + y + z - 1) - 2),
                    4*x*(3*x - y*(8*z - 3) + 3*y - 8*z*(x + y + z - 1) - 2),
                    32*x*y*(-x - y - 2*z + 1),
                ]
            )
        elif i == 5:  # between (1,2)
            phi = -4*y*(x + y + z - 1)*(8*x*z - 3*x - 3*z + 1)
            dphi = np.array(
                [
                    4*y*(-8*x*z + 3*x + 3*z - (8*z - 3)*(x + y + z - 1) - 1),
                    4*(-x - 2*y - z + 1)*(8*x*z - 3*x - 3*z + 1),
                    4*y*(-8*x*z + 3*x + 3*z - (8*x - 3)*(x + y + z - 1) - 1),
                ]
            )
        elif i == 6:  # between (0,2)
            phi = -4*x*(x + y + z - 1)*(8*y*z - 3*y - 3*z + 1)
            dphi = np.array(
                [
                    4*(-2*x - y - z + 1)*(8*y*z - 3*y - 3*z + 1),
                    4*x*(-8*y*z + 3*y + 3*z - (8*z - 3)*(x + y + z - 1) - 1),
                    4*x*(-8*y*z + 3*y + 3*z - (8*y - 3)*(x + y + z - 1) - 1),
                ]
            )
        elif i == 7:  # between (0,3)
            phi = 4*x*z*(3*x - 8*y*(x + y + z - 1) + 3*z - 2)
            dphi = np.array(
                [
                    4*z*(-x*(8*y - 3) + 3*x - 8*y*(x + y + z - 1) + 3*z - 2),
                    32*x*z*(-x - 2*y - z + 1),
                    4*x*(3*x - 8*y*(x + y + z - 1) - z*(8*y - 3) + 3*z - 2),
                ]
            )
        elif i == 8:
            phi = -4*y*z*(8*x*(x + y + z - 1) - 3*y - 3*z + 2)
            dphi = np.array(
                [
                    32*y*z*(-2*x - y - z + 1),
                    4*z*(-8*x*(x + y + z - 1) - y*(8*x - 3) + 3*y + 3*z - 2),
                    4*y*(-8*x*(x + y + z - 1) + 3*y - z*(8*x - 3) + 3*z - 2),
                ]
            )
        elif i == 9:
            phi = -4*z*(x + y + z - 1)*(8*x*y - 3*x - 3*y + 1)
            dphi = np.array(
                [
                    4*z*(-8*x*y + 3*x + 3*y - (8*y - 3)*(x + y + z - 1) - 1),
                    4*z*(-8*x*y + 3*x + 3*y - (8*x - 3)*(x + y + z - 1) - 1),
                    4*(-x - y - 2*z + 1)*(8*x*y - 3*x - 3*y + 1),
                ]
            )
        elif i == 10:
            phi = 27*x*y*(4*z - 1)*(x + y + z - 1)
            dphi = np.array(
                [
                   27*y*(4*z - 1)*(2*x + y + z - 1),
                   27*x*(4*z - 1)*(x + 2*y + z - 1),
                   27*x*y*(4*x + 4*y + 8*z - 5)
                ]
            )
        elif i == 11:
            phi = 27*x*y*z*(4*x + 4*y + 4*z - 3)
            dphi = np.array(
                [
                    27*y*z*(8*x + 4*y + 4*z - 3),
                    27*x*z*(4*x + 8*y + 4*z - 3),
                    27*x*y*(4*x + 4*y + 8*z - 3)
                ]
            )
        elif i == 12:
            phi = 27*y*z*(4*x - 1)*(x + y + z - 1)
            dphi = np.array(
                [
                  27*y*z*(8*x + 4*y + 4*z - 5),
                  27*z*(4*x - 1)*(x + 2*y + z - 1),
                  27*y*(4*x - 1)*(x + y + 2*z - 1)
                ]
            )
        elif i == 13:
            phi = 27*x*z*(4*y - 1)*(x + y + z - 1)
            dphi = np.array(
                [
                    27*z*(4*y - 1)*(2*x + y + z - 1),
                    27*x*z*(4*x + 8*y + 4*z - 5),
                    27*x*(4*y - 1)*(x + y + 2*z - 1)
                ]
            )
        elif i == 14:
            phi = 256*x*y*z*(-x - y - z + 1)
            dphi = np.array(
                [
                   256*y*z*(-2*x - y - z + 1),
                   256*x*z*(-x - 2*y - z + 1),
                   256*x*y*(-x - y - 2*z + 1)                
                ]
            )
        else:
            self._index_error()

        return phi, dphi

2D

class ElementTriCCR(ElementH1):

    nodal_dofs = 1
    facet_dofs = 1
    interior_dofs = 1
    maxdeg = 3
    dofnames = ['u', 'u', 'u']
    doflocs = np.array([[1., 0.],
                        [0., 1.],
                        [0., 0.],
                        [0.5, 0.5],
                        [0, 0.5],
                        [0.5, 0.],
                        [1./3, 1./3]])
    refdom = RefTri

    def lbasis(self, X, i):
        x, y = X

        if i == 0: 
            phi = -x*(-2*x + 3*y*(x + y - 1) + 1)
            dphi = np.array([6*x*y + 4*x + 3*y**2 - 3*y - 1,
                             3*x*(x + 2*y - 1)])
        elif i == 1:
            phi = -y*(3*x*(x + y - 1) - 2*y + 1)
            dphi = np.array([3*y*(-2*x - y + 1),
                             -3*x**2 - 6*x*y + 3*x + 4*y - 1])
        elif i == 2:
            phi = -(-2*x + y*(3*x - 2) + 1)*(x + y - 1)
            dphi = np.array([-6*x*y + 4*x - 3*y**2 + 7*y - 3,
                             -3*x**2 - 6*x*y + 7*x + 4*y - 3])
        elif i == 3:  # 0->1
            phi = 4*x*y*(3*x + 3*y - 2)
            dphi = np.array([4*y*(6*x + 3*y - 2),
                             4*x*(3*x + 6*y - 2)])
        elif i == 4:  # 1->2
            phi = 4*y*(3*x - 1)*(x + y - 1)
            dphi = np.array([4*y*(6*x + 3*y - 4),
                             4*(3*x - 1)*(x + 2*y - 1)])
        elif i == 5:  # 0->2
            phi = 4*x*(3*y - 1)*(x + y - 1)
            dphi = np.array([4*(3*y - 1)*(2*x + y - 1),
                             4*x*(3*x + 6*y - 4)])
        elif i == 6:
            phi = 27*x*y*(-x - y + 1)
            dphi = np.array(
                [
                    27*y*(-2*x - y + 1),
                    27*x*(-x - 2*y + 1)
                ]
            )
        else:
            self._index_error()

        return phi, dphi

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 have ElementTetDG and ElementTriDG already in scikit-fem, it seems that example 36 should be reproducible by simply changing the relevant spaces to

uelem = ElementVectorH1(ElementTetCCR())
pelem = ElementTetDG(ElementTetP1())

Please let me know what you think. In the meanwhile, I will check if I have implemented everything correctly (basis functions and their gradients)

bhaveshshrimali commented 3 years ago

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)

kinnala commented 3 years ago

Here are some tests that are run against the elements (if applicable):

Test Kronecker delta property of the nodes in doflocs: https://github.com/kinnala/scikit-fem/blob/master/tests/test_elements.py#L23
Test that derivatives give correct values: https://github.com/kinnala/scikit-fem/blob/master/tests/test_elements.py#L159
Test that the basis functions form a partition of unity: https://github.com/kinnala/scikit-fem/blob/master/tests/test_elements.py#L222

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: skfem/element/element_tri/element_tri_ccr.py and imports added to skfem/element/element_tri/__init__.py and skfem/element/__init__.py so that wildcard works properly. Then the above tests are modified accordingly.

If you want you can open a pull request with these changes.

kinnala commented 3 years ago

I can also do it if you want, it will be in July after I finish some other work.

bhaveshshrimali commented 3 years ago

Ok. I will test everything locally and then open a PR. We can take it from there. Thanks

bhaveshshrimali commented 3 years ago

Ok some progress...

Kronecker Delta

In [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 values

In [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 Unity

In [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 2D case. Will now get to checking example_36 with this to see if everything is indeed working as expected. Thanks

bhaveshshrimali commented 3 years ago

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...

bhaveshshrimali commented 3 years ago

I am able to reproduce example 36 now.

The mistake was in the dof ordering above. Changed it to match that in skfem.refdom and now everything seems to work. Will open a PR sometime over the weekend.

kinnala commented 3 years ago

Fantastic! I'll review the PR as soon as it's ready.

kinnala commented 3 years ago

Merged

kinnala / scikit-fem

Implementing Conforming CR #437 #651

3D

2D

Kronecker Delta

Derivative values

Partition of Unity