Formulators

Formulators specify the formulations of equation terms. The FormulationUP_pointwise stores the tensors for contraction, usually at integration points, with shape functions of the variables U and P. The Variables_<RDIM> contains the different geometrical variables of the context as well as the unknows variables U and P.

This a example of separation of concerns: the model and its application in a given context. The general form is:

template <> FormulationUP_pointwise apply_formulation(const MyPhysicalModel &, const Variables_<RDIM>& variables);

The model has generally some coefficients: for example the viscosity \(\nu\) in Stokes’s diffusion term. They are given under the form of expressions that may depend upon space or time. The formulator will translate them in real tensors (ScalarN) by contextualization coming from the variables, before building the formulation tensors.

Let’see first the composition of the Formulation_pointwise. Point-wise mean tensors are formed at the points where variables are defined: it is not yet the elementary vector and matrices, which will be later obtained automatically by the « applicationUP » framework. The goal is for the user to define only its custom part: the physical model action.

struct FormulationUP_pointwise
{
    uint dim;

    ScalarN   form0;  // energy/potential if it exists

    // For the elementary vector
    VectorN   formU0; // vector to be contracted with test functions of U
    MatrixN   formU1; // matrix to be contracted with test functions derivatives of U
    ScalarN   formP0; // scalar to be contracted with test functions of P
    VectorN   formP1; // vector to be contracted with test functions derivatives of P

    // For the elementary matrix: same but with a bilinear contraction, involving left and right test functions
    MatrixN   formU0_U0;
    O3MatrixN formU1_U0,formU0_U1;
    O4MatrixN formU1_U1;

    VectorN   formU0_P0,formP0_U0;
    MatrixN   formU1_P0,formU0_P1;
    O3MatrixN formU1_P1;

    ScalarN   formP0_P0;
    MatrixN   formP1_P1;
};

Here are the classical formulations for UP system, that will be used later.

// 1/2 |u|^2
template <> FormulationUP_pointwise apply_formulation < U,U > (const Variables_<RDIM>& variables)
{
    const uint dim = variables.dim;

    auto u = variables(_U_);

    return
    {
        .dim = dim,
        .form0     = 0.5*dot(u,u),
        .formU0    = u,
        .formU0_U0 = Tensor::identity_matrix(dim)
    };
}
// 1/2 |grad u|^2
template <> FormulationUP_pointwise apply_formulation < Grad<U>,Grad<U> > (const Variables_<RDIM>& variables)
{
    const uint dim = variables.dim;

    auto DuDx = variables.D(_U_);

    return
    {
        .dim = dim,
        .form0     = 0.5*ddot(DuDx,DuDx),
        .formU1    = DuDx,
        .formU1_U1 = Tensor::_I_(dim)
    };
}
// p div u (arising in the lagrangian)
template <> FormulationUP_pointwise apply_formulation < Div<U>,P > (const Variables_<RDIM>& variables)
{
    const uint dim = variables.dim;

    auto p      = variables  (_P_);
    auto DuDx   = variables.D(_U_);
    auto  div_u = trace(DuDx);
    auto I      = Tensor::identity_matrix(dim);

    return
    {
        .dim = dim,
        .form0     = p*div_u,
        .formU1    = p*I,
        .formP0    = div_u,
        .formU1_P0 = I
    };
}
// 1/2 |grad p|^2
template <> FormulationUP_pointwise apply_formulation < Grad<P>,Grad<P> > (const Variables_<RDIM>& variables)
{
    const uint dim = variables.dim;

    auto DpDx = variables.D(_P_);
    auto I    = Tensor::identity_matrix(dim);

    return
    {
        .dim = dim,
        .form0     = 0.5*dot(DpDx,DpDx),
        .formP1    = DpDx,
        .formP1_P1 = I
    };
}

Now physically-named formulators will be defined using the previous mathematically-named ones. They use models containing expressions for coefficients that will be traduced in real tensors (as some fields or functions) when applied in a context defining variables.

CFD terms

First of all, we have some physical models:

struct PressureStabilization : public OnDimension<RDIM>
{
    PressureStabilization(const Expr& b) : beta(b) {}
    Expr beta;
};
struct Poisson : public OnDimension<RDIM>
{
    Poisson(const Expr& n) : nu(n) {}
    Expr nu;
};
struct PoissonD : public OnDimension<RDIM>
{
    PoissonD(const Expr& n) : nu(n) {}
    Expr nu;
};
struct Stokes : public OnDimension<RDIM>
{
    Stokes(const Expr& n) : nu(n) {}
    Expr nu;
};
struct VelocityConvection : public OnDimension<RDIM>
{
};
struct NavierStokes : public OnDimension<RDIM>
{
    NavierStokes(const Expr& n, const Expr& b={}) : nu(n), beta(b) {}
    Expr nu;
    Expr beta;
};

For CFD, we can define readily:

template <> FormulationUP_pointwise apply_formulation(const Poisson& coefs, const Variables_<RDIM>& variables)
{
    auto nu = variables.scalarN(coefs.nu);

    return nu*apply_formulation<Grad<U>,Grad<U>> (variables);
}

template <> FormulationUP_pointwise apply_formulation(const Stokes& coefs, const Variables_<RDIM>& variables)
{
    auto nu = variables.scalarN(coefs.nu);
}
// the u.grad_u formulator
template <> FormulationUP_pointwise apply_formulation(const VelocityConvection&, const Variables_<RDIM>& variables)
{
    auto I = Tensor::identity_matrix(variables.dim);
    auto u    = variables  (_U_);
    auto DuDx = variables.D(_U_);
    return
    {
        .dim = variables.dim,
        .formU0    = dot(DuDx,u),
        .formU0_U0 = DuDx,
        .formU0_U1 = ...
    };
}
template <> FormulationUP_pointwise apply_formulation(const PressureStabilization& coefs, const Variables_<RDIM>& variables)
{
    if (!coefs.beta) return FormulationUP_pointwise{variables.dim};

    auto beta = variables.scalarN(coefs.beta);

    auto d = beta*variables.element_measure();

    return d*apply_formulation<Grad<P>,Grad<P>> (variables);
}
template <> FormulationUP_pointwise apply_formulation(const NavierStokes& coefs, const Variables_<RDIM>& variables)
{
    return apply_formulation(Stokes{coefs.nu},                 variables)
    +      apply_formulation(VelocityConvection{},             variables)
    +      apply_formulation(PressureStabilization{coefs.beta},variables);
}

Structural mechanics terms

First of all, the models:

struct ExternalPressure : public OnDimension<RCODIM1>
{
    ExternalPressure(const Expr& p) : external_pressure(p) {}
    Expr external_pressure;
};
template <RelativeDimension D> struct ExternalForce : public OnDimension<D>
{
    ExternalForce(const Expr& f) : external_force(f) {}
    Expr external_force;
};
template <RelativeDimension D> struct PenalizedPenetration : public OnDimension<D>
{
    PenalizedPenetration(const Expr& n, const Expr& l, const Expr& coef) : outward_normal(n/norm2(n)), level(l), K(coef) {}
    Expr outward_normal;
    Expr level;
    Expr K;
};

then the formulators:

template <> FormulationUP_pointwise apply_formulation(const Hyperelasticity& elasticity, const Variables_<RDIM>& variables)
{
    const uint dim = variables.dim;

    // all expressions of coefficients are replaced with tensors (ScalarN) using 'variables'
    auto material = apply(elasticity.material,variables);

    auto F = variables.F;

    if (elasticity.linearized)
    {
        auto Id = Tensor::identity_matrix(dim);
        auto Gu = F-Id;
        auto C  = Id+Gu+transpose(Gu);

        auto [W,WC,WCC] = Material::energyC_tensors(C,material,elasticity.kI);

        auto formulation = FormulationUP_pointwise
        {
            .dim       = dim,
            .form0     =   W,
            .formU1    = 2*WC,
            .formU1_U1 = 4*WCC
        };
        if (material.is_incompressible)
        {
            auto p = variables.get<0>(_P_);
            formulation +=
            {
                .form0      = p*trace(Gu),
                .formU1     = p*Id,
                .formP0     =   trace(Gu),
                .formU1_P0  =   Id
            };
        }
        return formulation;
    }

    auto [W,WF,WFF] = Material::energyF_tensors(F,material,elasticity.kI);

    assert(transpose(WFF)==WFF);

    auto formulation = FormulationUP_pointwise
    {
        .dim       = dim,
        .form0     = W,
        .formU1    = WF,
        .formU1_U1 = WFF
    };

    if (material.is_incompressible)
    {
        auto p = variables.get<0>(_P_);
        auto detF1 = det(F)-1;
        auto cofF  = cof(F);
        formulation +=
        {
            .form0     = p*detF1,
            .formU1    = p*cofF,
            .formU1_U1 = p*d2det(F),
            .formP0    = detF1,
            .formU1_P0 = cofF
        };
    }
    return formulation;
}

template <RelativeDimension D> FormulationUP_pointwise __apply(const ExternalForce<D>& spec, const Variables_<D>& variables)
{
    auto external_force = variables.vectorN(spec.external_force);

    return
    {
        .form0  = dot(external_force,variables.xF),
        .formU0 = external_force
    };
}
template <> FormulationUP_pointwise apply_formulation(const ExternalForce<RDIM>   & f, const Variables_<RDIM>   & v) { return __apply(f,v);}
template <> FormulationUP_pointwise apply_formulation(const ExternalForce<RDIM1>  & f, const Variables_<RDIM1>  & v) { return __apply(f,v);}
template <> FormulationUP_pointwise apply_formulation(const ExternalForce<RCODIM1>& f, const Variables_<RCODIM1>& v) { return __apply(f,v);}

template <RelativeDimension D> FormulationUP_pointwise apply_formulation(const PenalizedPenetration<D>& spec, const Variables_<D>& variables)
{
    auto x = variables.x;
    auto u = variables.u;
    auto xu = x+u;
    Decimal1 zero = 0;

    auto normal = variables.vector(spec.outward_normal);
    auto level0 = variables.scalar(spec.level);
    auto K      = variables.scalar(spec.K);

    auto level = dot(normal,xu)-level0;
    auto has_penetrated = test_lt0(level);

    auto Kn = select(has_penetrated,zero,K);

    return
    {
        .form0     = 0.5*K*pow(level,2)),
        .formU0    =     K*level            * normal,
        .formU0_U0 =     K * tensorial(normal,normal)
    };
}
template <> FormulationUP_pointwise apply_formulation(const ExternalPressure& spec, const Variables_<RCODIM1>& variables)
{
    auto external_pressure = variables.scalarN(spec.external_pressure);
    ...
}