Emily Boudreaux
cb20faeb3d
The test suite is extensive. Note that LLMs were used to aid in the development of the test suite. All LLM generated code was manually reviewed.
200 lines
8.1 KiB
C++
200 lines
8.1 KiB
C++
// Ported from MeanField tests/quadrature/policy.cpp using the consolidated form/role API.
|
|
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <format>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
|
|
#include <catch2/catch_test_macros.hpp>
|
|
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
|
#include <mfem.hpp>
|
|
|
|
#include "serif/backend/mfem/discretization/quadrature/mfem_resolver.hpp"
|
|
|
|
#include "serif/tests/test_tags.hpp"
|
|
|
|
namespace quadrature = serif::discretization::quadrature;
|
|
namespace forms = serif::discretization::forms;
|
|
namespace blocks = serif::discretization::blocks;
|
|
using VectorMass = forms::mass<blocks::displacement::vector, blocks::displacement::vector>;
|
|
using DensityMass = forms::mass<blocks::density::mass, blocks::density::mass>;
|
|
|
|
TEST_CASE(
|
|
"Quadrature Policy Computes Base Orders",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
const quadrature::Policy policy{quadrature::RuleSet{}};
|
|
quadrature::Query query{.form = forms::form_descriptor<VectorMass>};
|
|
query.trial_order = 3;
|
|
query.test_order = 3;
|
|
query.coefficient_order = 2;
|
|
query.geometry_weight_order = 1;
|
|
CHECK(policy.resolve(query).base_order == 9);
|
|
query.trial_order_reduction = 1;
|
|
query.test_order_reduction = 2;
|
|
CHECK(policy.resolve(query).base_order == 6);
|
|
query.trial_order_reduction = 5;
|
|
CHECK(policy.resolve(query).base_order == 4);
|
|
query.base_order = 7;
|
|
CHECK(policy.resolve(query).base_order == 7);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Quadrature Policy Composes Global Role And Form Boosts",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
quadrature::RuleSet rules;
|
|
rules.set_fallback_control({.boost = 1});
|
|
rules.set_role_control(quadrature::QuadratureRole::preconditioner, {.boost = 2});
|
|
rules.set_form_control<VectorMass>({.boost = 3});
|
|
const quadrature::Policy policy{rules};
|
|
const quadrature::Query query{
|
|
.form = forms::form_descriptor<VectorMass>, .role = quadrature::QuadratureRole::preconditioner, .base_order = 6
|
|
};
|
|
const auto result = policy.resolve(query);
|
|
CHECK(result.base_order == 6);
|
|
CHECK(result.boost == 6);
|
|
CHECK(result.order == 12);
|
|
CHECK_FALSE(result.used_fixed_order);
|
|
CHECK(policy.resolve({.form = forms::form_descriptor<DensityMass>, .base_order = 4}).order == 5);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Quadrature Fixed Orders Have Defined Precedence",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
quadrature::RuleSet rules;
|
|
const quadrature::Query query{.form = forms::form_descriptor<VectorMass>, .base_order = 6};
|
|
rules.set_fallback_control({.fixed_order = 9, .boost = 100});
|
|
CHECK(quadrature::Policy{rules}.resolve(query).order == 9);
|
|
rules.set_role_control(quadrature::QuadratureRole::discretization, {.fixed_order = 11, .boost = 200});
|
|
CHECK(quadrature::Policy{rules}.resolve(query).order == 11);
|
|
rules.set_form_control<VectorMass>({.fixed_order = 13, .boost = 300});
|
|
const auto result = quadrature::Policy{rules}.resolve(query);
|
|
CHECK(result.order == 13);
|
|
CHECK(result.used_fixed_order);
|
|
CHECK(result.boost == 0);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Quadrature Policy Rejects Invalid Orders",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
const quadrature::Policy policy{quadrature::RuleSet{}};
|
|
quadrature::Query query{.form = forms::form_descriptor<VectorMass>};
|
|
query.base_order = -1;
|
|
CHECK_THROWS_AS(policy.resolve(query), std::invalid_argument);
|
|
query.base_order.reset();
|
|
query.trial_order = -1;
|
|
CHECK_THROWS_AS(policy.resolve(query), std::invalid_argument);
|
|
query.trial_order = 1;
|
|
query.test_order_reduction = -1;
|
|
CHECK_THROWS_AS(policy.resolve(query), std::invalid_argument);
|
|
quadrature::RuleSet rules;
|
|
rules.set_form_control<VectorMass>({.fixed_order = -1});
|
|
query.test_order_reduction = 0;
|
|
CHECK_THROWS_AS(quadrature::Policy{rules}.resolve(query), std::invalid_argument);
|
|
rules.set_form_control<VectorMass>({.boost = -2});
|
|
CHECK_THROWS_AS(quadrature::Policy{rules}.resolve(query), std::invalid_argument);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"MFEM Rule Factory Includes Geometry And Dynamic Orders",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
using namespace quadrature::backend::mfem;
|
|
const MFEMResolver resolver{quadrature::Policy{quadrature::RuleSet{}}};
|
|
const std::array dynamic_orders{2, 1};
|
|
MFEMQuadratureContext context{
|
|
.geometry = mfem::Geometry::CUBE, .geometry_weight_order = 1, .dynamic_orders = dynamic_orders
|
|
};
|
|
const auto rule = resolver.resolve<VectorMass>(context);
|
|
CHECK(rule.resolution.base_order == 10);
|
|
CHECK(rule.integration_rule == &mfem::IntRules.Get(mfem::Geometry::CUBE, 10));
|
|
context.geometry_weight_order = -1;
|
|
CHECK_THROWS_AS(resolver.resolve<VectorMass>(context), std::invalid_argument);
|
|
context.geometry_weight_order = 0;
|
|
const std::array invalid_orders{-1};
|
|
context.dynamic_orders = invalid_orders;
|
|
CHECK_THROWS_AS(resolver.resolve<VectorMass>(context), std::invalid_argument);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"MFEM Quadrature Rule Integrates Tensor Polynomial Exactly",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
const quadrature::backend::mfem::MFEMResolver resolver{quadrature::Policy{quadrature::RuleSet{}}};
|
|
const auto rule = resolver.resolve<VectorMass>({.geometry = mfem::Geometry::CUBE});
|
|
double integral = 0.0;
|
|
for (int i = 0; i < rule.integration_rule->GetNPoints(); ++i) {
|
|
const auto& point = rule.integration_rule->IntPoint(i);
|
|
integral += point.weight * std::pow(point.x, 4) * std::pow(point.y, 2) * point.z;
|
|
}
|
|
CHECK_THAT(integral, Catch::Matchers::WithinAbs(1.0 / 30.0, 1.0e-14));
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Policy Controlled Vector Mass Assembly Matches Overintegrated Reference",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
auto mesh = mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 1.2, 0.9, 1.4);
|
|
mfem::H1_FECollection collection(3, 3);
|
|
mfem::FiniteElementSpace space(&mesh, &collection, 3);
|
|
auto& transformation = *mesh.GetElementTransformation(0);
|
|
const auto& element = *space.GetFE(0);
|
|
const quadrature::backend::mfem::MFEMResolver resolver{quadrature::Policy{quadrature::RuleSet{}}};
|
|
mfem::VectorMassIntegrator controlled;
|
|
const auto resolution = resolver.configure<VectorMass>(
|
|
controlled, {.geometry = element.GetGeomType(), .geometry_weight_order = transformation.OrderW()}
|
|
);
|
|
mfem::VectorMassIntegrator reference;
|
|
reference.SetIntRule(&mfem::IntRules.Get(element.GetGeomType(), resolution.order + 8));
|
|
mfem::DenseMatrix actual, expected;
|
|
controlled.AssembleElementMatrix(element, transformation, actual);
|
|
reference.AssembleElementMatrix(element, transformation, expected);
|
|
actual -= expected;
|
|
CHECK(actual.FNorm() / expected.FNorm() < 1.0e-12);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Integration Rule Adapter Supports Each Advertised Signature",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
struct ReferenceIntegrator {
|
|
const mfem::IntegrationRule* selected = nullptr;
|
|
void SetIntegrationRule(const mfem::IntegrationRule& rule) {
|
|
selected = &rule;
|
|
}
|
|
} reference;
|
|
struct PointerIntegrator {
|
|
const mfem::IntegrationRule* selected = nullptr;
|
|
void SetIntRule(const mfem::IntegrationRule* rule) {
|
|
selected = rule;
|
|
}
|
|
} pointer;
|
|
struct IntRuleReferenceIntegrator {
|
|
const mfem::IntegrationRule* selected = nullptr;
|
|
void SetIntRule(const mfem::IntegrationRule& rule) {
|
|
selected = &rule;
|
|
}
|
|
} int_reference;
|
|
|
|
const quadrature::backend::mfem::MFEMResolver resolver{quadrature::Policy{quadrature::RuleSet{}}};
|
|
const quadrature::backend::mfem::MFEMQuadratureContext context{.geometry = mfem::Geometry::CUBE};
|
|
resolver.configure<VectorMass>(reference, context);
|
|
resolver.configure<VectorMass>(pointer, context);
|
|
resolver.configure<VectorMass>(int_reference, context);
|
|
CHECK(reference.selected == &mfem::IntRules.Get(mfem::Geometry::CUBE, 6));
|
|
CHECK(pointer.selected == reference.selected);
|
|
CHECK(int_reference.selected == reference.selected);
|
|
}
|
|
|
|
TEST_CASE(
|
|
"Discrete Functional Space Can Be Formatted",
|
|
tags::unit & tags::quadrature
|
|
) {
|
|
const blocks::displacement::vector::FunctionalSpace space;
|
|
CHECK(std::format("{}", space).find("rank=Vector") != std::string::npos);
|
|
}
|