Flexible integration utility

[btalamini]

Introduce a new utility to compute integrals and evaluate functions at quadrature points that provides a bunch of features missing from the FunctionSpace class. The main goal is flexibility:

• Integrate/evaluate functions with an arbitrary number of input fields
• Integrals/evaluations differentiable on all input fields, including the mesh coordinates. That is, do shape optimization without constructing a new FunctionSpace object each design iteration.
• Input fields of different types:
  • standard finite element nodal fields
  • DG nodal fields
  • quadrature fields (e.g. internal variables)
  • Something I'm calling "uniform fields", things that have a single value for the whole mesh (like time)
• Ability to specify the integrand signature, in terms of field values and field gradients. That is, handle different cases like $f(u)$, or $f(\nabla u, u)$, or $f(u, \nabla u, v, t)$
• Extensibility, so we can add new field types with different interpolation and gradient operators

It's easiest to explain with some examples.

Let's implement the energy for the Helmholtz equation:

$$I(u) = \int_{B} \frac{1}{2} \left( u^2 + \nabla u \cdot \nabla u \right) \, dv$$

where $u$ is a field in the standard finite element space of piecewise second order polynomials.

# Choose function spaces for the inputs
# Here we pick a standard Lagrange polynomial space of degree 2
input_spaces = PkField(2, self.mesh),
# The energy density function has the form f(u, grad(u))
FIELD_INDEX = 0
integrand_signature = Value(FIELD_INDEX), Gradient(FIELD_INDEX)
# Instantiate the entity to compute the integral
field_operator = FieldOperator(spaces, integrand_signature, mesh, quad_rule)

# Define the function to integrate.
# It must have the signature corresponding to integrand_signature
def f(u, dudX):
    return 0.5*(u*u + np.dot(dudX, dxdX))

energy = field_opearator.integrate(f, mesh.coords, U)

Another example: linear elasticity, with a DG $P_0$ shear modulus field:

# Choose function spaces for the inputs
input_spaces = PkField(2, self.mesh), DG_PkField(0, self.mesh)
DISPLACEMENT = 0
SHEAR_MODULUS = 1
integrand_signature = Gradient(DISPLACEMENT), Value(SHEAR_MODULUS)
field_operator = FieldOperator(spaces, integrand_signature, mesh, quad_rule)

def f(dudX, mu):
    lam = 1.0
    strain = sym(dudX)
    return mu*np.tensordot(strain, strain) + lam*np.trace(strain)**2

# Code would be here to create the input fields with specific values,
# U = (something)
# mu = (something else)

energy = field_opearator.integrate(f, mesh.coords, U, mu)

To-do in following PRs:

• Extend to surface integrals/evaluations
• Use FieldOperator to re-write the Mechanics and other physics modules.

[codecov]

Codecov Report

❌ Patch coverage is 88.64865% with 21 lines in your changes missing coverage. Please review.
✅ Project coverage is 76.06%. Comparing base (14311a7) to head (acc95ff).
⚠️ Report is 40 commits behind head on main.

Files with missing lines	Patch %	Lines
optimism/FieldOperator.py	88.64%	21 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #120      +/-   ##
==========================================
+ Coverage   75.56%   76.06%   +0.50%     
==========================================
  Files          65       66       +1     
  Lines        5647     5832     +185     
==========================================
+ Hits         4267     4436     +169     
- Misses       1380     1396      +16     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[cmhamel]

@ralberd docs removed from CI