warping function as poisson's equation

[numenic]

Hello there,

I try to deep in scikit-fem. I have trouble to solve what should be an easy problem. I'm dealing with structural beams cross-sections, eg. shear center calculation.

Manuals say that solving the omega function over the domain is the way to go. Here is a derivation (from https://engineering.stackexchange.com/a/6721):

.

The equation is a simple laplacian, therefore:

    @skfem.BilinearForm
    def laplace(u, v, w):
        return dot(grad(u), grad(v))

For the second equation, well... here is my attempt, but it is obviously wrong:

    @skfem.LinearForm
    def L(v, w):
        x, y = w.x  # global coordinates
        f = (x**2 + y**2) * 0.5
        return f * v

I would really appreciate some help here...

[kinnala]

Because of the tangential derivative, I think you either want to integrate-by-parts or then you can try to interpolate x^2 + y^2 at integration points and write the tangential derivative with the help of a normal vector. I can write you an example later today.

[kinnala]

I tried replicating one of the results from you link.

Here is the code. I'm using the version from master, hopefully it also works in your version.

from skfem import *
from skfem.models import laplace
from skfem.helpers import dot, grad
import numpy as np

m1 = MeshQuad.init_tensor(
    np.linspace(0, 200, 41),
    np.linspace(0, 10, 3),
)

m2 = MeshQuad.init_tensor(
    np.linspace(0, 200, 41),
    np.linspace(290, 300, 3),
)

m3 = MeshQuad.init_tensor(
    np.linspace(95, 105, 3),
    np.linspace(10, 290, 40),
)


m = m1 + m2 + m3
# it seems you need to move the origin
# so that it matches the center of mass
# to get similar result as at stackexchange
m = m.translated((-100, -150))

basis = Basis(m, ElementQuad2())
fbasis = basis.boundary()

A = laplace.assemble(basis)

@LinearForm
def linf(v, w):
    s = np.array([-w.n[1], w.n[0]])
    return dot(grad(w['yz']), s)

X = basis.global_coordinates()
yz = basis.project(.5 * (X[0] ** 2 + X[1] ** 2))

b = linf.assemble(fbasis, yz=fbasis.interpolate(yz))

u = solve(*condense(A, b, D=m.nodes_satisfying(lambda x: (x[0]==100)*(x[1]==150))))

basis.plot(u, shading='gouraud', colorbar=True).show()

[kinnala]

Perhaps you want to change the solve line to

u = solve(*condense(A, b, D=m.nodes_satisfying(lambda x: (x[0]==0)*(x[1]==0))))

so that you are fixing the same node to zero as in the reference.

[numenic]

wow. Thanks. I try to understand what you did. What are the coordinates stored in X ? values do not look like "real" coordinates, although the number of coordinates match the number of elements...

[kinnala]

X contains the locations of integration points.