Need more documentation about TractionType

[tepperly]

https://github.com/LLNL/serac/blob/eb736b5abe45304aacc987ce6004286f721ab692/src/serac/physics/solid_mechanics.hpp#L496

I need more documentation about how the TractionType lambda needs to be defined (standard disclaimers about not being a mechanical engineer). I kind of expected that x would be a ::mfem::Vector indicating the location of the point in space, and I expected that its size would be equal to the spatial dimension of the problem. Based on difficult to interpret compiler output, it appears to be a ::serac::tensor<double,2> which leaves me scratching my head about what x is here.

The second argument is named normal, so I would expect that it's a ::mfem::Vector of size equal to the spatial dimension. Presumably, it's pointing in the direction from inside to outside? I am not currently using normal, so I don't know if these interpretations are wrong or right. Ideally, the documentation would make it clearer.

Our LiDO test that compiles, has a return value of 1*x, so presumably, the return value is supposed to be a ::serac::tensor<double, 2>. If x is the spatial location in some sense, why would the load be in the direction of the spatial location?

In the example I am working on, I am trying to set a point external load in the negative y direction at the top in y and in the middle in x for a three point bend problem. My incorrect intuitions led me to right this (in the LiDO API call that calls setPiolaTraction):

    physics_->setBoundaryLoad([X_WIDTH,Y_WIDTH,X_LEN,Y_LEN](const auto &location, [[maybe_unused]] const auto &normal)
    {
        constexpr double LOAD=-0.5;
        ::mfem::Vector result(2);
        const double signedDistanceFromXCenter(location[0]-X_LEN*0.5);
        const double signedDistanceFromYTop(Y_LEN-location[1]);
        result[0] = 0.0; // never any x-component                                                                                                                                      
        result[1] = ((::std::fabs(signedDistanceFromXCenter) < X_WIDTH) &&
                     (::std::fabs(signedDistanceFromYTop) < Y_WIDTH)) ? LOAD : 0.0;
        return result;
    });

There also appears to be a cut-and-paste error when the current documentation says "return the thermal flux value".

[tepperly]

Is this a correct way to do what I was hoping?

[X_WIDTH,Y_WIDTH,X_LEN,Y_LEN](const auto &location, [[maybe_unused]] const auto &normal)
    {
        if constexpr (::std::is_same<decltype(location), const ::serac::tensor<double,2> &>::value)
        {
            constexpr double LOAD=-0.5;
            const double signedDistanceFromXCenter(location[0]-X_LEN*0.5);
            const double signedDistanceFromYTop(Y_LEN-location[1]);
            return (((::std::fabs(signedDistanceFromXCenter) < X_WIDTH) &&
                     (::std::fabs(signedDistanceFromYTop) < Y_WIDTH)) ? LOAD : 0.0)*normal;
        }
        else
        {
            // derivative is zero                                                                                                                                                      
            return 0*normal;
        }
    }

[jamiebramwell]

This (and the material class interface) definitely needs documentation. An example of this type of class is shown here:

https://github.com/LLNL/serac/blob/d5dcd67d7bb717b5f8337766a7854f2c058240b1/src/serac/physics/materials/solid_material.hpp#L202-L220

We don't use mfem::Vector for these low level types and instead use serac::tensor for compatibility with AD dual types. Given your code snippet, I think this is what you want:

[X_WIDTH,Y_WIDTH,X_LEN,Y_LEN](const auto &location, [[maybe_unused]] const auto &normal, double t)
    {
            // this is needed for shape-aware physics
            auto traction = 0.0 * location + 0.0 * normal;
            constexpr double LOAD=-0.5;
            const double signedDistanceFromXCenter(location[0]-X_LEN*0.5);
            const double signedDistanceFromYTop(Y_LEN-location[1]);

            if (((::std::fabs(signedDistanceFromXCenter) < X_WIDTH) && (::std::fabs(signedDistanceFromYTop) < Y_WIDTH)) {
                  traction += LOAD * normal;
            }

            return traction;
    }

[samuelpmishLLNL]

I need more documentation about how the TractionType lambda needs to be defined (standard disclaimers about not being a mechanical engineer).

I'll write this here first, to answer your immediate question (and we'll follow up soon, improving the documentation with more information about the types passed in to q-functions). As for function signatures, there are broadly two kinds:

• "q-functions" for domain integrals get passed (in order):

  1. the spatial position of the quadrature point (e.g. tensor<double, dim>) as the first argument
  2. a variadic list of tuples, containing the {value, derivative} of each trial space field evaluated at the quadrature point

• "q-functions" for boundary integrals get passed:

  1. the spatial position of the quadrature point (e.g. tensor<double, dim>) as the first argument
  2. the outward-facing unit normal of the surface, at the quadrature point (e.g. tensor<double, dim>) as the second argument
  3. a variadic list of tuples, containing the {value, derivative} of each trial space field evaluated at the quadrature point

The time value at the end of these kernels is a confusing temporary inclusion. That part of the interface is going to change soon (to better support differentiating with respect to time and other quantities that aren't finite element fields).

User: So why do the examples use auto for all the q-function arguments?

This is a language limitation of C++17-- lambda function arguments have to be concrete types (e.g. tensor<double, dim>) or auto. Unfortunately, concrete types do not admit differentiation, so we have to use auto, which makes it less clear to the user what types are actually being passed.

I kind of expected that x would be a ::mfem::Vector indicating the location of the point in space

That's the correct expectation mathematically, but mfem::Vector is not the right container to use in a kernel for a variety of reasons:

• its size is unknown
• it allocates dynamically memory
• it does not support automatic differentiation
• it does not work on a GPU
• it's missing a lot of useful functions (e.g. operator overloads, etc)

For reasons like this, we pass tensors of the appropriate shape to these q-functions.

serac::tensor<double,2> which leaves me scratching my head about what x is here.

It's exactly what you expected-- serac::tensor<double,2> is a vector of 2 components

Our LiDO test that compiles, has a return value of 1*x, so presumably, the return value is supposed to be a ::serac::tensor<double, 2>.

Boundary q-functions only return the "source" terms, as described here. For H1 and L2 spaces, the source terms are vectors with as many components as the test space was defined with (e.g. normal heat flux for thermal, tractions for elasticity).

In the example I am working on, I am trying to set a point external load in the negative y direction at the top in y and in the middle in x for a three point bend problem. My incorrect intuitions led me to right this (in the LiDO API call that calls setPiolaTraction): ...

I'd recommend that LIDO (or maybe serac) define some additional types for doing basic geometric queries like:

// put these definitions somewhere

// axis-aligned bounding box
template < int dim >
struct AABB {
  tensor< double, dim > min, max;
};

// is the point x inside the given bounding box?
template < typename T, int dim >
bool inside(tensor<T, dim> x, AABB<dim> bounds){
  for (int i = 0; i < dim; i++) {
    if (x[i] < bounds.min[i] || x[i] > bounds.max[i]) return false;
  }
  return true;
}

Then, your q-function is much simpler to write and understand:

AABB<2> bounds {...};
[bounds](const auto & x, const auto & normal, double t) {
  constexpr double LOAD=-0.5;
  return inside(x, bounds) * LOAD * normal;
}

Note: this sort of geometric query is not differentiable w.r.t to spatial position

---

edit: make inside() implementation support different types stored in the position vector

[tepperly]

Is ::serac::tensor<type, dim> always an order one tensor (i.e., something that one might call a vector)?

[samuelpmishLLNL]

Is ::serac::tensor<type, dim> always an order one tensor (i.e., something that one might call a vector)?

Yes, see: https://serac.readthedocs.io/en/latest/sphinx/dev_guide/tensor_dual.html

tensor< T, n > is a vector of n components tensor< T, n1, n2 > is an n1-by-n2, matrix of Ts tensor< T, n1, n2, n3 > represents a n1-by-n2-by-n3 (and so on) tensor of Ts and so on

[tepperly]

Based on compiler output, it appears that x can be either a ::serac::tensor<double, 2> or a serac::tensor<serac::dual<serac::tuple<serac::tensor<double, 2>, serac::zero>>, 2>. It doesn't always have a operator[] method defined, so I don't think the approach with inside will work.