how to customize the class ElasticityIntegrator

[xrfe]

Hi mfem developer, I want to modify the ElasticityIntegrator to incorperate the thermal stress [−α(3λ+2μ)(T−T0​)I], below is my mofification, however, I got no convergency. Do you know what goes wrong? Thanks!

class ThermoElasticityIntegrator : public ElasticityIntegrator
{
protected:
    Coefficient α  // Coefficient of thermal expansion
    Coefficient &temperature;  // Temperature field as a Coefficient
    double T0;  // Reference temperature

public:
    ThermoElasticityIntegrator(Coefficient &lambda, Coefficient &mu, 
                               Coefficient &alpha_, Coefficient &temp, 
                               double T0_)
        : ElasticityIntegrator(lambda, mu), alpha(alpha_), 
          temperature(temp), T0(T0_) { }

    virtual void AssembleElementMatrix(const FiniteElement &el,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat) override
    {
        int dof = el.GetDof(), dim = el.GetDim();
      
        DenseMatrix B(dim*dim, dof*dim);
        DenseMatrix DB(dim*dim, dof*dim);
        DenseMatrix D(dim*dim, dim*dim);
        Vector shape(dof);
        DenseMatrix dshape(dof, dim);
        DenseMatrix gshape(dof, dim);

        elmat.SetSize(dof*dim);
        elmat = 0.0;

        const IntegrationRule *ir = IntRule;
        if (ir == NULL)
        {
            ir = &IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 1);
        }

        for (int i = 0; i < ir->GetNPoints(); i++)
        {
            const IntegrationPoint &ip = ir->IntPoint(i);
            Trans.SetIntPoint(&ip);
            el.CalcDShape(ip, dshape);
            el.CalcShape(ip, shape);

            MultAtB(dshape, Trans.InverseJacobian(), gshape);
         
            double w = ip.weight * Trans.Weight();
            double lambda_val = lambda->Eval(Trans, ip);
            double mu_val = mu->Eval(Trans, ip);
            double alpha_val = alpha.Eval(Trans, ip);

            // Compute thermal strain
            double T = temperature.Eval(Trans, ip);  // Use Eval method for Coefficient
            double thermal_strain = alpha_val * (T - T0);

            // Assemble B matrix
            B = 0.0;
            for (int j = 0; j < dof; j++)
            {
                for (int k = 0; k < dim; k++)
                {
                    for (int l = 0; l < dim; l++)
                    {
                        B(k*dim+l, j*dim+k) = gshape(j, l);
                    }
                }
            }

            // Compute D matrix (plane strain case for 2D, or 3D)
            D = 0.0;
            double lambda_2mu = lambda_val + 2.0 * mu_val;
            for (int k = 0; k < dim; k++)
            {
                D(k*dim+k, k*dim+k) = lambda_2mu;
                for (int l = 0; l < dim; l++)
                {
                    if (k != l)
                    {
                        D(k*dim+k, l*dim+l) = lambda_val;
                        D(k*dim+l, k*dim+l) = mu_val;
                    }
                }
            }

            // Apply thermal strain
            Vector thermal_strain_vec(dim*dim);
            for (int j = 0; j < dim*dim; j++)
            {
                thermal_strain_vec(j) = thermal_strain;
            }

            Vector thermal_stress(dim*dim);
            D.Mult(thermal_strain_vec, thermal_stress);

            // Compute DB = D * B
            Mult(D, B, DB);

            // Assemble element matrix
            AddMult_a_ABt(w, B, DB, elmat);

            // Add thermal stress contribution
            for (int j = 0; j < dof; j++)
            {
                for (int k = 0; k < dim; k++)
                {
                    for (int l = 0; l < dim; l++)
                    {
                        elmat(j*dim+k, j*dim+l) -= w * shape(j) * thermal_stress(k*dim+l);
                    }
                }
            }
        }
    }
};

[v-dobrev]

It seems like the thermal stress does not depend on the deformation, so it should not be added to the ElasticityIntegrator. Instead it defines a known force that can be added to the r.h.s. of the equation. The linear form should be $(\nabla\cdot \sigma_t, \vec{v})$ where $\sigma_t$ is the thermal stress and $\vec{v}$ is the test function. Since $\sigma_t = \beta I$ the linear form can be written as $(\nabla \beta, \vec{v})$. Integrating by parts this becomes:
$$(\nabla \beta, \vec{v}) = \int_{\Omega} \nabla \beta \cdot \vec{v} dx = -\int_{\Omega} \beta \nabla\cdot \vec{v} dx + \int_{\Gamma_N} \beta \vec{v}\cdot\vec{n} dS$$
For the second (boundary) term one can use the integrator VectorBoundaryFluxLFIntegrator. For the first term, however, it looks like we don't have implementation for this specific linear form (with test functions in vector $H^1$-conforming FE space), so it will need to be implemented. We can help with that, if the above seems reasonable approach for you to use.

[xrfe]

Thanks for the reponses. That's a very valuable suggestion.
A quicker resolution would be:

1. calculate the gradient of the temperature field,
2. use the computed gradient as a force term in the elasticity problem with VectorDomainLFIntegrator.

[v-dobrev]

Yes, if $\beta$ has no jumps (e.g. due to jumps in material properties) then computing its gradient and using it as a force should work.

If you represent $\beta$ as an $H^1$-conforming (Par)GridFunction, you can use the GradientGridFunctionCoefficient class to perform the gradient computation.

[j-signorelli]

Here is code for a custom LinearFormIntegrator I have written and verified in development of my own linear elasticity solver with thermal expansion, $(f, \nabla\cdot\vec{v})$ for $H^1$:

.hpp:

class VectorDomainLFDivIntegrator : public mfem::LinearFormIntegrator
{   
    protected:
        mfem::Coefficient &fcoeff;
        mfem::DenseMatrix dshape, gshape;
        mfem::Vector pelvect;
    public:
        VectorDomainLFDivIntegrator(mfem::Coefficient &fcoeff_) : LinearFormIntegrator(), fcoeff(fcoeff_) {};
        void AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &Tr, mfem::Vector &elvect);
};

.cpp:

void VectorDomainLFDivIntegrator::AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
{
    int dof = el.GetDof();
    int dim = el.GetDim();
    double c, w;

    dshape.SetSize(dof, dim);
    gshape.SetSize(dof, dim);

    elvect.SetSize(dof * dim);
    pelvect.SetSize(dof * dim);

    const IntegrationRule *ir = IntRule;
    if (ir == NULL)
    {
        int order = 2 * Tr.OrderGrad(&el);
        ir = &IntRules.Get(el.GetGeomType(), order);
    }
    
    elvect = 0.0;
    pelvect = 0.0;

    for (int i = 0; i < ir -> GetNPoints(); i++)
    {
        const IntegrationPoint &ip = ir->IntPoint(i);
    
        el.CalcDShape(ip, dshape);
    
        Tr.SetIntPoint(&ip);
        c = fcoeff.Eval(Tr,ip);
        w = c*ip.weight*Tr.Weight();
        Mult(dshape, Tr.InverseJacobian(), gshape);
        gshape.GradToDiv(pelvect);

        pelvect *= w;
        elvect += pelvect;
    }
}

@v-dobrev would there be interest/value in submitting a PR for this VectorDomainLFDivIntegrator ?

[xrfe]

Thanks for sharing the codes, yes it works and generates identical results as the method suggested by @v-dobrev.
It's also great if mfem can extend the VectorFEDomainLFDivIntegrator to support more spaces.
BTW, I got warning when comipling it using nvcc, not sure if you experience the following message:
warning #611-D: overloaded virtual function "mfem::LinearFormIntegrator::AssembleRHSElementVect" is only partially overridden in class "VectorDomainLFDivIntegrator"