Components of the boundary stress

[gittsoy]

I have a question related to issue #618. I'm trying to compute a normal component of the boundary stress vector:

We can compute components s[itr, jtr] and project to ElementQuad1() as suggested in #618. Then for the rectangular area, the value of the normal component will coincide with orthogonal normal stress. The question is what to do when the area is not a rectangle.

I tried to use the following code:

sigma = linear_stress(Lambda, mu)
u_bottom = bottom_basis.interpolate(u)
eps_gp = 0.5 * (u_bottom.grad + u_bottom.grad.transpose(1, 0, 2, 3))
sig_gp = sigma(eps_gp)
normals=bottom_basis.normals.value
sigma_n = (sig_gp[0][0]*normals[0]+sig_gp[0][1]*normals[1])*normals[0] \ 
         +(sig_gp[1][0]*normals[0]+sig_gp[1][1]*normals[1])*normals[1]

I obtain 2D numpy.ndarray with size [number of facets at the border x 5] (I'm using ElementQuad2() element for displacements).
How can i obtain 1D array with the values at the boundary dofs? In other words i need some function reverse for 'interpolate' or smth like this.

Sorry for bothering and thanks for previous answers

[kinnala]

M = BilinearForm(lambda u, v, w: u * v).assemble(bottom_basis)
f = LinearForm(lambda v, w: sigma_n * v).assemble(bottom_basis)
sigma_n_dofs = solve(M, f)

Can you try it?

[gittsoy]

I tried it but got an error while doing assemble:

~/.conda/envs/fem/lib/python3.8/site-packages/skfem/assembly/form/bilinear_form.py in assemble(self, ubasis, vbasis, **kwargs)
     71                 rows[ixs] = vbasis.element_dofs[i]
     72                 cols[ixs] = ubasis.element_dofs[j]
---> 73                 data[ixs] = self._kernel(
     74                     ubasis.basis[j],
     75                     vbasis.basis[i],

ValueError: could not broadcast input array from shape (2,5) into shape (120,)

Here is code for mesh and basis

m = MeshQuad.init_tensor(np.linspace(0, 3, 121), np.linspace(0, 1, 41))
e1 = ElementQuad2()
e = ElementVectorH1(e1)
gb = InteriorBasis(m, e, intorder=4)

bottom_facets = m.facets_satisfying(lambda x: x[1] == 0.)
bottom_basis = FacetBasis(m, e, facets=bottom_facets)

[kinnala]

Alright, I was a bit confused. Maybe do

bottom_scalar_basis = FacetBasis(m, e1, facets=bottom_facets)

and then as above:

M = BilinearForm(lambda u, v, w: u * v).assemble(bottom_scalar_basis)
f = LinearForm(lambda v, w: sigma_n * v).assemble(bottom_scalar_basis)
sigma_n_dofs = solve(M, f)

I don't have the full code so cannot test it out. Anyhow, the idea is to perform a projection onto the finite element basis.

[gittsoy]

It didn't work. I've got an error: "linsolve.py:206: MatrixRankWarning: Matrix is exactly singular warn("Matrix is exactly singular", MatrixRankWarning)". Btw i'm using 2.4 version from conda. Here is full code that reproduces the problem:

import numpy as np
from skfem import *
from skfem.models.elasticity import *
from skfem.visuals.matplotlib import *
from skfem.helpers import dot

def tracRight(x, y):
    pr = 1.
    return np.array([0., pr])

@LinearForm
def loadingN_Right(v, w):
    return dot(tracRight(*w.x), v)

m = MeshQuad.init_tensor(np.linspace(0, 3, 121), np.linspace(0, 1, 41))
e1 = ElementQuad2()
e = ElementVectorH1(e1)
gb = InteriorBasis(m, e, intorder=4)

bottom_facets = m.facets_satisfying(lambda x: x[1] == 0.)
bottom_basis = FacetBasis(m, e, facets=bottom_facets)
right_facets = m.facets_satisfying(lambda x: x[0] == 3.)
right_basis = FacetBasis(m, e, facets=right_facets)

clamped = gb.get_dofs(lambda x: x[0] == 0.0).all()
free = gb.complement_dofs(clamped)
loaded = gb.get_dofs(lambda x: x[1] == 2.0).all()

Lambda, mu = lame_parameters(1000.0, 0.3)
K = asm(linear_elasticity(Lambda, mu), gb)
rpN = asm(loadingN_Right, right_basis)
u = solve(*condense(K, rpN, I=free)) 

sigma = linear_stress(Lambda, mu)
u_bottom = bottom_basis.interpolate(u)
eps_gp = 0.5 * (u_bottom.grad + u_bottom.grad.transpose(1, 0, 2, 3))
sig_gp = sigma(eps_gp)
normals=bottom_basis.normals.value
sigma_n = (sig_gp[0][0]*normals[0]+sig_gp[0][1]*normals[1])*normals[0]+(sig_gp[1][0]*normals[0]+sig_gp[1][1]*normals[1])*normals[1]

bottom_scalar_basis = FacetBasis(m, e1, facets=bottom_facets)
M = BilinearForm(lambda u, v, w: u * v).assemble(bottom_scalar_basis)
f = LinearForm(lambda v, w: sigma_n * v).assemble(bottom_scalar_basis)
sigma_n_dofs = solve(M, f)

As result, i got array sigma_n_dofs with size 19521.

[kinnala]

Alright, I forgot that M has all the columns and rows for composability with other discrete operators and is not invertible unless condense is used. I'll try to fix it.

[kinnala]

This is the modification I made:

sigma_n_dofs = solve(*condense(M, f, I=bottom_scalar_basis.get_dofs(bottom_facets)))

Now sigma_n_dofs has one value per P2 node in the original mesh. The dofs not belonging to bottom_facets should be zero. If you want to use P1 nodes instead, you must use ElementQuad1 when defining bottom_scalar_basis and make sure that it has the same quadrature points as bottom_basis.

[gittsoy]

Thanks, think i got it. As i understood i just need to pick values from sigma_n_dofs corresponding to boundary dofs. Here is code i used:

import matplotlib
import matplotlib.pyplot as plt

bottom_dofs = bottom_scalar_basis.get_dofs(lambda x: x[1] == 0.).all()
gb_scalar = InteriorBasis(m, e1, intorder=4)

dofs = gb_scalar.get_dofs(bottom_facets)
x = gb_scalar.doflocs[:, dofs.flatten()][0]
normal_stress = sigma_n_dofs[bottom_dofs]
p = x.argsort()
normal_stress = normal_stress[p]
x = x[p]

plt.rcParams["figure.figsize"] = (7,3)
plt.plot(x,normal_stress)

[gittsoy]

If you want to use P1 nodes instead, you must use ElementQuad1 when defining bottom_scalar_basis and make sure that it has the same quadrature points as bottom_basis.

Thus do i only need to change the definition of bottom_scalar_basis using ElementQuad1 and not bottom_basis itself? I tried:

bottom_scalar_basis = FacetBasis(m, ElementQuad1(), intorder=4, facets=bottom_facets)

but got an error:

<ipython-input-28-da4fa9d94ce1> in <lambda>(v, w)
      8 bottom_scalar_basis = FacetBasis(m, ElementQuad1(), intorder=4, facets=bottom_facets)
      9 M = BilinearForm(lambda u, v, w: u * v).assemble(bottom_scalar_basis)
---> 10 f = LinearForm(lambda v, w: sigma_n * v).assemble(bottom_scalar_basis)
     11 
     12 sigma_n_dofs = solve(*condense(M, f, I=bottom_scalar_basis.get_dofs(bottom_facets)))

ValueError: operands could not be broadcast together with shapes (120,5) (120,3) 

How can i make them have the same quadrature points? Or i misunderstood and just need to change e1 = ElementQuad1()?
If i do this it's working and i obtain the more smooth graph of normal stress.

[kinnala]

You can init both bases with same intorder.

[gittsoy]

Thanks for the help, I finally figured it out.