Coefficient defined on quadrature points in the form definition

[tianyikillua]

I have a the coefficient function that does not belong to the same finite element space as v, but is defined for each quadrature point.

@LinearForm
def loading(v, w):
    return dot(w["sig0"], sym(grad(v)))

Is it possible?

Background
In linear elasticity, sometimes the stress - strain is given by

sigma = Hooke @ eps + sig0

where sig0 is the initial stress tensor, possibly heterogeneous (different value for each Gauss point). This can be used for instance to model thermal stress.

[tianyikillua]

Oh, I see DiscreteField meets exactly what I asked. Is there an example which illustrates how numerically this field can be defined, i.e., its shape and ordering?

[gdmcbain]

A DiscreteField is usually generated by InteriorBasis.interpolate

scikit-fem/skfem/assembly/basis/basis.py

Lines 182 to 183 in 72ba058

	def interpolate(self, w: ndarray) -> Union[DiscreteField,
	Tuple[DiscreteField, ...]]:

I think a similar issue to the coefficient function not belonging to the same finite element space arises in mixed finite element methods as used for Stokes equations (exx 18, 24, 27, 30)

scikit-fem/docs/examples/ex18.py

Lines 56 to 59 in 91a276a

	element = {'u': ElementVectorH1(ElementTriP2()),
	'p': ElementTriP1()}
	basis = {variable: InteriorBasis(mesh, e, intorder=3)
	for variable, e in element.items()}

I think it doesn't matter that that the spaces are different so long as the quadrature points coincide, i.e. the InteriorBasis that you interpolate from doesn't have to be the one of your unknown.

Can you create an InteriorBasis for your initial stress, making sure to use the same quadrature points, populate it, e.g. using skfem.utils.project, and then interpolate?

One recently introduced way to make a new InteriorBasis on the same Mesh with the same quadrature rule is the InteriorBasis.with_element method.

[tianyikillua]

Thanks for your reply!

I see. This means simply if ElementTriP1 is used for my displacement, I only need to define my sig0 as a ElementTriP0 field, right?

Now the problem is that my sig0 is not a scalar but a (symmetric) d-dim tensor, to be double contracted by sym(grad(v))). Do you have a workround for this?

[gdmcbain]

Yes, I thought of this after writing and have been thinking about it. I've only tried what I wrote above for scalars and vectors, not a tensor like stress; there isn't a higher order equivalent of ElementVector. I wonder whether one's required for something like initial stress. It might depend on where it comes from (a previous boundary value problem?). If it's thermal stress, so that it's derived from a scalar temperature field, then I'm thinking not, it should be possible to proceed as in the machinery in skfem.models.elasticity in which for example a vectorial displacement field is turned into a second-order with sym_grad.

If the thermal stress is an isotropic tensor proportional to the excess temperature, it might be appropriate to use skfem.helpers.eye, as used in the λ part of the elastic linear_stress.

[tianyikillua]

Absolutely...Unfortunately my initial stress is neither expressible by sym_grad, nor is a diagonal tensor. It's coming from another physics which I calculate elsewhere...I can only provide the numerical values of sig0 (3 components for 2-d elasticity).

Do you think it is possible to tweak ElementVector by regarding self.dim as the number of components, and not the spatial dimension?

[kinnala]

I think no tweaks are required in this case:

import numpy as np
from skfem import *
from skfem.helpers import sym_grad, ddot

mesh = MeshTri()
basis = InteriorBasis(mesh, ElementVector(ElementTriP1()))

Z = basis.zero_w()  # has shape Nelems x Nqp
sig0 = DiscreteField(np.array([[Z, Z], [Z, Z]]))  # change Z to something nonzero

@LinearForm
def loading(v, w):
    return ddot(w["sig0"], sym_grad(v))
    
b = loading.assemble(basis, sig0=sig0)

In long term I have a plan to include something like ElementTensor. If you want to try implementing it, it should be pretty similar to how ElementVector is done.

[kinnala]

Let me add that if you need the global coordinates corresponding to each element of Z, you can get them as follows:

In [6]: x, y = basis.global_coordinates().value                                 

In [7]: x                                                                       
Out[7]: 
array([[0.16666667, 0.66666667, 0.16666667],
       [0.83333333, 0.33333333, 0.83333333]])

In [8]: y                                                                       
Out[8]: 
array([[0.16666667, 0.16666667, 0.66666667],
       [0.33333333, 0.83333333, 0.83333333]])

It's not exactly clear to me how you are planning to calculate the components of sig0 so there might also be a simpler way.

Edit: Also make sure to fix the version of scikit-fem because these hints are relying on somewhat low level details that might change more often.

[tianyikillua]

I'm testing the "data-driven computational mechanics" approach, where I need to prescribe initial stress fields computed by a local (at every Gauss point) spatial tree query in the phase space (eps, sig).

Thanks for the code snippet, so I suppose I could do something like

basis = InteriorBasis(mesh, ElementTriP0()))  # I think `ElementVector` is not needed right?

sig0_11 = basis.zero_w()  # has shape Nelems x Nqp
sig0_11[:, :] = ...
# same for 22 and 12
sig0 = DiscreteField(np.array([[sig0_11, sig0_12], [sig0_12, sig0_22]])

for my displacement in ElementTriP1.

[kinnala]

Actually, the element used in InteriorBasis depends on what you want to v be. If you specify it as ElementTriP0() then sym_grad(v) does not make sense I think? It has nothing to do with sig0 other than that you get the correct shape from basis.zero_w().

Edit: Maybe this part of the tutorial will help you or maybe not: https://scikit-fem.readthedocs.io/en/latest/forms.html#forms-return-numpy-arrays

Edit 2: If you use two different InteriorBasis objects then you must manually verify that they use the same integration points, e.g., by explicitly specifying the integration order: InteriorBasis(..., intorder=4)

[tianyikillua]

Oh I understand better now! Thank you. I will report if I have something.

[tianyikillua]

Sorry I may have a quick but related question. From a solution vector v, is there is post-processing procedure with which I can evaluate its gradient (the stress / strain fields) on Gauss points?

[tianyikillua]

I found this code snippet from ex04.py, and it's working properly I think.

element_sig = ElementTriDG(ElementTriP1())
basis_sig = InteriorBasis(mesh, element_sig)

@BilinearForm
def proj_stress(u, v, w):
    return sigma(sym_grad(u))[0, 0] * v

@BilinearForm
def mass(u, v, w):
    return u * v

solve(asm(mass, basis_sig), asm(proj_strain, basis, basis_sig) @ u)  # for sig_11

I wonder if there is a more "local" way to do this without doing a global solve (like LocalSolver in FEniCS)...

[kinnala]

I wonder if there is a more "local" way to do this without doing a global solve (like LocalSolver in FEniCS)...

Unfortunately, no. I think the "correct" way to do this would be to return here a list of local stiffness matrices:

scikit-fem/skfem/assembly/form/bilinear_form.py

Line 101 in 91a276a

# TODO: allow user to change, e.g. cuda or petsc

And then use np.linalg.solve or similar which supports solving multiple dense systems at one go. I can try implementing it in few days if you want.

[tianyikillua]

Yes I think it could be useful. However note that here we need in fact the local B(n_elem, n_gp_elem, n_strains, n_dof_local) matrices that compute the local strain tensors (or for Poisson, the gradient) for each Gauss points...with this we should be able to parallelize the computation of local gradients during post-processing.

[kinnala]

I wasn't thinking anything complicated; in scikit-fem we aim to focus on the assembly part and leave parallelization to other tools. I just wanted to point out that a call to skfem.utils.solve could be replaced by a call to np.linalg.solve (which is probably parallel if there are multiple small dense matrices) because asm(mass, basis_sig) is block-diagonal and we should not need to assemble those into one large matrix at all.

[gdmcbain]

If you pass v to the .interpolate method of its basis, you'll get a DiscreteField which has a .grad attribute.

[tianyikillua]

It appears to be working fine. Thanks all!

[QsingularityAi]

Is there any way to write 1D and 2D for nonlinear elastic constitutive relations using skfem, I am totally stuck. how to define bilinear form for this.