492 lines
25 KiB
C++
492 lines
25 KiB
C++
|
#include <algorithm>
|
||
|
|
#include <cmath>
|
||
|
|
#include <limits>
|
||
|
|
#include <stdexcept>
|
||
|
|
#include <utility>
|
||
|
|
|
||
|
|
#include <catch2/catch_test_macros.hpp>
|
||
|
|
|
||
|
|
#include <mfem.hpp>
|
||
|
|
#include <mpi.h>
|
||
|
|
|
||
|
|
import mean_field;
|
||
|
|
import test_helpers;
|
||
|
|
|
||
|
|
namespace radial_extension_test_utils {
|
||
|
|
namespace deformation = mean_field::deformation;
|
||
|
|
namespace domain = mean_field::utils::domain;
|
||
|
|
namespace field = mean_field::field;
|
||
|
|
|
||
|
|
using Schema = domain::CoreEnvelopeVacuumDomainSchema;
|
||
|
|
|
||
|
|
[[nodiscard]] mfem::Vector referenceCenter(const int spatialDimension) {
|
||
|
|
mfem::Vector center(spatialDimension);
|
||
|
|
center = 0.0;
|
||
|
|
return center;
|
||
|
|
}
|
||
|
|
|
||
|
|
[[nodiscard]] deformation::PreparedNodalRadialSurface makePreparedSurface(const mean_field::fem::FEM &fem) {
|
||
|
|
const field::ScalarBoundaryDofMap surfaceDofMap =
|
||
|
|
field::make_stellar_surface_scalar_dof_map<Schema>(*fem.surfaceDeformationFes);
|
||
|
|
const deformation::SurfaceDeformationCompilationContext context{*fem.surfaceDeformationFes, surfaceDofMap};
|
||
|
|
return deformation::compileSurfaceDeformationPrescription(
|
||
|
|
deformation::NodalRadialSurface{referenceCenter(fem.mesh->SpaceDimension())}, context
|
||
|
|
);
|
||
|
|
}
|
||
|
|
|
||
|
|
[[nodiscard]] double globalInnerProduct(
|
||
|
|
const mfem::Vector &first,
|
||
|
|
const mfem::Vector &second,
|
||
|
|
MPI_Comm communicator
|
||
|
|
) {
|
||
|
|
REQUIRE(first.Size() == second.Size());
|
||
|
|
const double local = first * second;
|
||
|
|
double global = 0.0;
|
||
|
|
MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, communicator);
|
||
|
|
return global;
|
||
|
|
}
|
||
|
|
|
||
|
|
[[nodiscard]] double relativeError(
|
||
|
|
const mfem::Vector &actual,
|
||
|
|
const mfem::Vector &expected
|
||
|
|
) {
|
||
|
|
REQUIRE(actual.Size() == expected.Size());
|
||
|
|
mfem::Vector difference(actual);
|
||
|
|
difference -= expected;
|
||
|
|
return difference.Norml2() / std::max(expected.Norml2(), std::numeric_limits<double>::epsilon());
|
||
|
|
}
|
||
|
|
|
||
|
|
[[nodiscard]] int mfemByNodesVectorDof(
|
||
|
|
const int scalarTrueDof,
|
||
|
|
const int component,
|
||
|
|
const int scalarTrueDofCount
|
||
|
|
) {
|
||
|
|
return scalarTrueDof + component * scalarTrueDofCount;
|
||
|
|
}
|
||
|
|
} // namespace radial_extension_test_utils
|
||
|
|
|
||
|
|
TEST_CASE(
|
||
|
|
"Radial Interior And Vacuum Extensions Advertise Closed Form Boundary Behavior",
|
||
|
|
tags::radial_deformation_extension_validation
|
||
|
|
) {
|
||
|
|
namespace deformation = mean_field::deformation;
|
||
|
|
|
||
|
|
STATIC_CHECK(deformation::InteriorDeformationExtension<deformation::PowerLawRadialInteriorExtension>);
|
||
|
|
STATIC_CHECK(
|
||
|
|
deformation::PreparedInteriorDeformationExtension<deformation::PreparedPowerLawRadialInteriorExtension>
|
||
|
|
);
|
||
|
|
STATIC_CHECK(deformation::VacuumDeformationExtension<deformation::FixedInfinityRadialVacuumExtension>);
|
||
|
|
STATIC_CHECK(
|
||
|
|
deformation::PreparedVacuumDeformationExtension<deformation::PreparedFixedInfinityRadialVacuumExtension>
|
||
|
|
);
|
||
|
|
|
||
|
|
const deformation::PowerLawRadialInteriorExtension interior;
|
||
|
|
const deformation::InteriorDeformationExtensionDescriptor interiorDescriptor = interior.descriptor();
|
||
|
|
CHECK(interior.radialPower() == 2.0);
|
||
|
|
CHECK(interiorDescriptor.name == "PowerLawRadialInteriorExtension");
|
||
|
|
CHECK(interiorDescriptor.linearOnReferenceGeometry);
|
||
|
|
CHECK(interiorDescriptor.requiresRadialFoliation);
|
||
|
|
CHECK_FALSE(interiorDescriptor.requiresAuxiliarySolve);
|
||
|
|
CHECK(interiorDescriptor.supportsExactNewtonLinearization());
|
||
|
|
CHECK(interiorDescriptor.centerBehavior == deformation::InteriorCenterBehavior::FixedAtReferenceCenter);
|
||
|
|
|
||
|
|
const deformation::FixedInfinityRadialVacuumExtension vacuum;
|
||
|
|
const deformation::VacuumDeformationExtensionDescriptor vacuumDescriptor = vacuum.descriptor();
|
||
|
|
CHECK(vacuumDescriptor.name == "FixedInfinityRadialVacuumExtension");
|
||
|
|
CHECK(vacuumDescriptor.linearOnReferenceGeometry);
|
||
|
|
CHECK(vacuumDescriptor.requiresRadialFoliation);
|
||
|
|
CHECK_FALSE(vacuumDescriptor.requiresAuxiliarySolve);
|
||
|
|
CHECK(vacuumDescriptor.supportsExactNewtonLinearization());
|
||
|
|
CHECK(vacuumDescriptor.outerBoundaryBehavior == deformation::VacuumOuterBoundaryBehavior::FixedAtReferenceInfinity);
|
||
|
|
|
||
|
|
CHECK_THROWS_AS(deformation::PowerLawRadialInteriorExtension{0.5}, std::invalid_argument);
|
||
|
|
CHECK_THROWS_AS(
|
||
|
|
deformation::PowerLawRadialInteriorExtension{std::numeric_limits<double>::infinity()}, std::invalid_argument
|
||
|
|
);
|
||
|
|
}
|
||
|
|
|
||
|
|
TEST_CASE(
|
||
|
|
"Radial Extensions Reproduce The Stellar Surface Fix Reference Infinity And Preserve Positive Volume Maps",
|
||
|
|
tags::radial_deformation_extension_analytic &tags::radial_deformation_extension_mapping
|
||
|
|
) {
|
||
|
|
namespace deformation = mean_field::deformation;
|
||
|
|
namespace domain = mean_field::utils::domain;
|
||
|
|
namespace field = mean_field::field;
|
||
|
|
|
||
|
|
mean_field::utils::Args args = test_utils::setup_args();
|
||
|
|
mean_field::fem::FEM fem = mean_field::fem::setup_fem(args.mesh_file, args, 0);
|
||
|
|
REQUIRE(fem.okay());
|
||
|
|
|
||
|
|
const deformation::PreparedNodalRadialSurface surface = radial_extension_test_utils::makePreparedSurface(fem);
|
||
|
|
const deformation::RadialDeformationExtensionCompilationContext context =
|
||
|
|
deformation::makeRadialDeformationExtensionCompilationContext<radial_extension_test_utils::Schema>(
|
||
|
|
*fem.surfaceDeformationFes, *fem.displacementFes, *fem.logicalReferenceMesh
|
||
|
|
);
|
||
|
|
const deformation::PreparedPowerLawRadialInteriorExtension interior =
|
||
|
|
deformation::compileInteriorDeformationExtension(deformation::PowerLawRadialInteriorExtension{}, context);
|
||
|
|
const deformation::PreparedPowerLawRadialInteriorExtension cubicInterior =
|
||
|
|
deformation::compileInteriorDeformationExtension(deformation::PowerLawRadialInteriorExtension{3.0}, context);
|
||
|
|
const deformation::PreparedFixedInfinityRadialVacuumExtension vacuum =
|
||
|
|
deformation::compileVacuumDeformationExtension(deformation::FixedInfinityRadialVacuumExtension{}, context);
|
||
|
|
|
||
|
|
REQUIRE(interior.surfaceDisplacementSize() == surface.surfaceDisplacementSize());
|
||
|
|
REQUIRE(vacuum.surfaceDisplacementSize() == surface.surfaceDisplacementSize());
|
||
|
|
REQUIRE(interior.interiorDisplacementSize() == fem.displacementFes->GetTrueVSize());
|
||
|
|
REQUIRE(vacuum.vacuumDisplacementSize() == fem.displacementFes->GetTrueVSize());
|
||
|
|
|
||
|
|
constexpr double surfaceAmplitude = 0.02;
|
||
|
|
mfem::Vector parameters(surface.parameterCount());
|
||
|
|
parameters = surfaceAmplitude;
|
||
|
|
mfem::Vector surfaceDisplacement(surface.surfaceDisplacementSize());
|
||
|
|
surface.buildSurfaceDisplacement(parameters, surfaceDisplacement);
|
||
|
|
|
||
|
|
mfem::Vector interiorDisplacement(interior.interiorDisplacementSize());
|
||
|
|
mfem::Vector vacuumDisplacement(vacuum.vacuumDisplacementSize());
|
||
|
|
interior.buildInteriorDisplacement(surfaceDisplacement, interiorDisplacement);
|
||
|
|
vacuum.buildVacuumDisplacement(surfaceDisplacement, vacuumDisplacement);
|
||
|
|
|
||
|
|
const int spatialDimension = fem.mesh->SpaceDimension();
|
||
|
|
const field::ScalarBoundaryDofMap stellarSurfaceMap =
|
||
|
|
field::make_scalar_boundary_dof_map<domain::StellarSurface, radial_extension_test_utils::Schema>(
|
||
|
|
*fem.surfaceDeformationFes
|
||
|
|
);
|
||
|
|
const field::ScalarBoundaryDofMap infinitySurfaceMap =
|
||
|
|
field::make_scalar_boundary_dof_map<domain::InfinitySurface, radial_extension_test_utils::Schema>(
|
||
|
|
*fem.surfaceDeformationFes
|
||
|
|
);
|
||
|
|
|
||
|
|
CHECK_THROWS_AS(
|
||
|
|
deformation::RadialDeformationExtensionCompilationContext(
|
||
|
|
*fem.surfaceDeformationFes, *fem.displacementFes, *fem.logicalReferenceMesh, stellarSurfaceMap,
|
||
|
|
stellarSurfaceMap,
|
||
|
|
domain::make_attribute_marker<domain::Stellar, radial_extension_test_utils::Schema>(*fem.mesh),
|
||
|
|
domain::make_attribute_marker<domain::Vacuum, radial_extension_test_utils::Schema>(*fem.mesh),
|
||
|
|
radial_extension_test_utils::Schema::template boundary_attribute<domain::StellarSurface>(),
|
||
|
|
radial_extension_test_utils::Schema::template boundary_attribute<domain::StellarSurface>()
|
||
|
|
),
|
||
|
|
std::invalid_argument
|
||
|
|
);
|
||
|
|
|
||
|
|
CHECK_THROWS_AS(
|
||
|
|
deformation::makeRadialDeformationExtensionCompilationContext<radial_extension_test_utils::Schema>(
|
||
|
|
*fem.surfaceDeformationFes, *fem.displacementFes, *fem.mesh
|
||
|
|
),
|
||
|
|
std::invalid_argument
|
||
|
|
);
|
||
|
|
|
||
|
|
const double stellarSurfaceRadius = context.stellarSurfaceLogicalRadius();
|
||
|
|
const double infinitySurfaceRadius = context.infinitySurfaceLogicalRadius();
|
||
|
|
REQUIRE(stellarSurfaceRadius > 0.0);
|
||
|
|
REQUIRE(infinitySurfaceRadius > stellarSurfaceRadius);
|
||
|
|
|
||
|
|
constexpr double tolerance = 2.0e-11;
|
||
|
|
bool hasInterpolatedSurfacePoint = false;
|
||
|
|
for (int scalarDof = 0; scalarDof < interior.scalarTrueDofCount(); ++scalarDof) {
|
||
|
|
const double referenceRadius = context.logicalRadius(scalarDof);
|
||
|
|
const int interpolationEntryCount = context.surfaceInterpolationEntryCount(scalarDof);
|
||
|
|
|
||
|
|
if (interpolationEntryCount == 0) {
|
||
|
|
CHECK(referenceRadius <= 64.0 * std::numeric_limits<double>::epsilon() * infinitySurfaceRadius);
|
||
|
|
} else {
|
||
|
|
double interpolationWeightSum = 0.0;
|
||
|
|
for (int entry = 0; entry < interpolationEntryCount; ++entry) {
|
||
|
|
const int surfaceCoordinate = context.surfaceGlobalCoordinate(scalarDof, entry);
|
||
|
|
CHECK(surfaceCoordinate >= 0);
|
||
|
|
CHECK(surfaceCoordinate < stellarSurfaceMap.global_size());
|
||
|
|
interpolationWeightSum += context.surfaceInterpolationWeight(scalarDof, entry);
|
||
|
|
}
|
||
|
|
CHECK(std::abs(interpolationWeightSum - 1.0) <= tolerance);
|
||
|
|
hasInterpolatedSurfacePoint |= interpolationEntryCount > 1;
|
||
|
|
}
|
||
|
|
|
||
|
|
if (interior.hasStellarSupport(scalarDof)) {
|
||
|
|
const double expectedWeight =
|
||
|
|
referenceRadius == 0.0 ? 0.0 : std::pow(referenceRadius / stellarSurfaceRadius, 2.0);
|
||
|
|
CHECK(std::abs(interior.radialWeight(scalarDof) - expectedWeight) <= tolerance);
|
||
|
|
const double expectedCubicWeight =
|
||
|
|
referenceRadius == 0.0 ? 0.0 : std::pow(referenceRadius / stellarSurfaceRadius, 3.0);
|
||
|
|
CHECK(std::abs(cubicInterior.radialWeight(scalarDof) - expectedCubicWeight) <= tolerance);
|
||
|
|
} else {
|
||
|
|
CHECK(interior.radialWeight(scalarDof) == 0.0);
|
||
|
|
for (int component = 0; component < spatialDimension; ++component) {
|
||
|
|
const int volumeVectorDof = radial_extension_test_utils::mfemByNodesVectorDof(
|
||
|
|
scalarDof, component, interior.scalarTrueDofCount()
|
||
|
|
);
|
||
|
|
CHECK(interiorDisplacement(volumeVectorDof) == 0.0);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (vacuum.hasVacuumSupport(scalarDof)) {
|
||
|
|
const double expectedWeight =
|
||
|
|
(infinitySurfaceRadius - referenceRadius) / (infinitySurfaceRadius - stellarSurfaceRadius);
|
||
|
|
CHECK(std::abs(vacuum.radialWeight(scalarDof) - expectedWeight) <= tolerance);
|
||
|
|
} else {
|
||
|
|
CHECK(vacuum.radialWeight(scalarDof) == 0.0);
|
||
|
|
for (int component = 0; component < spatialDimension; ++component) {
|
||
|
|
const int volumeVectorDof = radial_extension_test_utils::mfemByNodesVectorDof(
|
||
|
|
scalarDof, component, vacuum.scalarTrueDofCount()
|
||
|
|
);
|
||
|
|
CHECK(vacuumDisplacement(volumeVectorDof) == 0.0);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
CHECK(hasInterpolatedSurfacePoint);
|
||
|
|
|
||
|
|
const double componentValues[3]{1.25, -0.75, 2.5};
|
||
|
|
mfem::Vector constantSurfaceDisplacement(surface.surfaceDisplacementSize());
|
||
|
|
for (int surfaceDof = 0; surfaceDof < surface.parameterCount(); ++surfaceDof) {
|
||
|
|
for (int component = 0; component < spatialDimension; ++component) {
|
||
|
|
constantSurfaceDisplacement(spatialDimension * surfaceDof + component) = componentValues[component];
|
||
|
|
}
|
||
|
|
}
|
||
|
|
mfem::Vector constantInteriorDisplacement(interior.interiorDisplacementSize());
|
||
|
|
interior.buildInteriorDisplacement(constantSurfaceDisplacement, constantInteriorDisplacement);
|
||
|
|
|
||
|
|
mfem::Vector radialWeightTrueDofs(interior.scalarTrueDofCount());
|
||
|
|
for (int scalarDof = 0; scalarDof < interior.scalarTrueDofCount(); ++scalarDof) {
|
||
|
|
radialWeightTrueDofs(scalarDof) =
|
||
|
|
interior.hasStellarSupport(scalarDof) ? interior.radialWeight(scalarDof) : 0.0;
|
||
|
|
}
|
||
|
|
mfem::ParGridFunction radialWeightField(fem.surfaceDeformationFes.get());
|
||
|
|
mfem::ParGridFunction constantVectorField(fem.displacementFes.get());
|
||
|
|
radialWeightField.SetFromTrueDofs(radialWeightTrueDofs);
|
||
|
|
constantVectorField.SetFromTrueDofs(constantInteriorDisplacement);
|
||
|
|
|
||
|
|
const mfem::Array<int> stellarMarker =
|
||
|
|
domain::make_attribute_marker<domain::Stellar, radial_extension_test_utils::Schema>(*fem.mesh);
|
||
|
|
int sampledStellarElement = -1;
|
||
|
|
for (int element = 0; element < fem.mesh->GetNE() && sampledStellarElement < 0; ++element) {
|
||
|
|
const int attribute = fem.mesh->GetAttribute(element);
|
||
|
|
if (attribute > 0 && attribute <= stellarMarker.Size() && stellarMarker[attribute - 1] != 0) {
|
||
|
|
sampledStellarElement = element;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
REQUIRE(sampledStellarElement >= 0);
|
||
|
|
const mfem::IntegrationPoint &samplePoint =
|
||
|
|
mfem::Geometries.GetCenter(fem.mesh->GetElementBaseGeometry(sampledStellarElement));
|
||
|
|
const double sampledRadialWeight = radialWeightField.GetValue(sampledStellarElement, samplePoint);
|
||
|
|
mfem::Vector sampledVector(spatialDimension);
|
||
|
|
constantVectorField.GetVectorValue(sampledStellarElement, samplePoint, sampledVector);
|
||
|
|
for (int component = 0; component < spatialDimension; ++component) {
|
||
|
|
CHECK(std::abs(sampledVector(component) - componentValues[component] * sampledRadialWeight) <= tolerance);
|
||
|
|
}
|
||
|
|
|
||
|
|
mfem::Vector arbitrarySurfaceDisplacement(surface.surfaceDisplacementSize());
|
||
|
|
for (int dof = 0; dof < arbitrarySurfaceDisplacement.Size(); ++dof) {
|
||
|
|
const double index = static_cast<double>(dof + 1);
|
||
|
|
arbitrarySurfaceDisplacement(dof) = 0.03 * std::sin(0.29 * index) - 0.01 * std::cos(0.17 * index);
|
||
|
|
}
|
||
|
|
mfem::Vector arbitraryInteriorDisplacement(interior.interiorDisplacementSize());
|
||
|
|
mfem::Vector arbitraryVacuumDisplacement(vacuum.vacuumDisplacementSize());
|
||
|
|
interior.buildInteriorDisplacement(arbitrarySurfaceDisplacement, arbitraryInteriorDisplacement);
|
||
|
|
vacuum.buildVacuumDisplacement(arbitrarySurfaceDisplacement, arbitraryVacuumDisplacement);
|
||
|
|
|
||
|
|
for (int surfaceDof = 0; surfaceDof < stellarSurfaceMap.local_size(); ++surfaceDof) {
|
||
|
|
const int scalarDof = stellarSurfaceMap.volume_true_dof(surfaceDof);
|
||
|
|
for (int component = 0; component < spatialDimension; ++component) {
|
||
|
|
const int surfaceVectorDof = spatialDimension * surfaceDof + component;
|
||
|
|
const int volumeVectorDof =
|
||
|
|
radial_extension_test_utils::mfemByNodesVectorDof(scalarDof, component, interior.scalarTrueDofCount());
|
||
|
|
CHECK(
|
||
|
|
std::abs(
|
||
|
|
arbitraryInteriorDisplacement(volumeVectorDof) - arbitrarySurfaceDisplacement(surfaceVectorDof)
|
||
|
|
) <= tolerance
|
||
|
|
);
|
||
|
|
CHECK(
|
||
|
|
std::abs(
|
||
|
|
arbitraryVacuumDisplacement(volumeVectorDof) - arbitrarySurfaceDisplacement(surfaceVectorDof)
|
||
|
|
) <= tolerance
|
||
|
|
);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
for (int infinityDof = 0; infinityDof < infinitySurfaceMap.local_size(); ++infinityDof) {
|
||
|
|
const int scalarDof = infinitySurfaceMap.volume_true_dof(infinityDof);
|
||
|
|
for (int component = 0; component < spatialDimension; ++component) {
|
||
|
|
const int volumeVectorDof =
|
||
|
|
radial_extension_test_utils::mfemByNodesVectorDof(scalarDof, component, vacuum.scalarTrueDofCount());
|
||
|
|
CHECK(std::abs(vacuumDisplacement(volumeVectorDof)) <= tolerance);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
mfem::Vector combinedDisplacement(interiorDisplacement);
|
||
|
|
for (int scalarDof = 0; scalarDof < vacuum.scalarTrueDofCount(); ++scalarDof) {
|
||
|
|
if (!vacuum.hasVacuumSupport(scalarDof) || interior.hasStellarSupport(scalarDof)) {
|
||
|
|
continue;
|
||
|
|
}
|
||
|
|
for (int component = 0; component < spatialDimension; ++component) {
|
||
|
|
const int volumeVectorDof =
|
||
|
|
radial_extension_test_utils::mfemByNodesVectorDof(scalarDof, component, vacuum.scalarTrueDofCount());
|
||
|
|
combinedDisplacement(volumeVectorDof) = vacuumDisplacement(volumeVectorDof);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
mfem::ParGridFunction displacement(fem.displacementFes.get());
|
||
|
|
displacement.SetFromTrueDofs(combinedDisplacement);
|
||
|
|
double localMinimumDeterminant = std::numeric_limits<double>::infinity();
|
||
|
|
for (int element = 0; element < fem.mesh->GetNE(); ++element) {
|
||
|
|
mfem::ElementTransformation *transformation = fem.mesh->GetElementTransformation(element);
|
||
|
|
const mfem::FiniteElement *finiteElement = fem.displacementFes->GetFE(element);
|
||
|
|
const mfem::IntegrationRule &rule =
|
||
|
|
mfem::IntRules.Get(transformation->GetGeometryType(), finiteElement->GetOrder() + 2);
|
||
|
|
|
||
|
|
for (int point = 0; point < rule.GetNPoints(); ++point) {
|
||
|
|
transformation->SetIntPoint(&rule.IntPoint(point));
|
||
|
|
mfem::DenseMatrix displacementGradient;
|
||
|
|
displacement.GetVectorGradient(*transformation, displacementGradient);
|
||
|
|
for (int component = 0; component < spatialDimension; ++component) {
|
||
|
|
displacementGradient(component, component) += 1.0;
|
||
|
|
}
|
||
|
|
localMinimumDeterminant = std::min(localMinimumDeterminant, displacementGradient.Det());
|
||
|
|
}
|
||
|
|
}
|
||
|
|
double globalMinimumDeterminant = 0.0;
|
||
|
|
MPI_Allreduce(&localMinimumDeterminant, &globalMinimumDeterminant, 1, MPI_DOUBLE, MPI_MIN, fem.mesh->GetComm());
|
||
|
|
CHECK(globalMinimumDeterminant > 0.0);
|
||
|
|
}
|
||
|
|
|
||
|
|
TEST_CASE(
|
||
|
|
"Radial Extension Jacobians Match Centered Differences And Their Transposes Preserve Virtual Work",
|
||
|
|
tags::radial_deformation_extension_linearization
|
||
|
|
) {
|
||
|
|
namespace deformation = mean_field::deformation;
|
||
|
|
|
||
|
|
mean_field::utils::Args args = test_utils::setup_args();
|
||
|
|
mean_field::fem::FEM fem = mean_field::fem::setup_fem(args.mesh_file, args, 0);
|
||
|
|
REQUIRE(fem.okay());
|
||
|
|
|
||
|
|
const deformation::RadialDeformationExtensionCompilationContext context =
|
||
|
|
deformation::makeRadialDeformationExtensionCompilationContext<radial_extension_test_utils::Schema>(
|
||
|
|
*fem.surfaceDeformationFes, *fem.displacementFes, *fem.logicalReferenceMesh
|
||
|
|
);
|
||
|
|
const deformation::PreparedPowerLawRadialInteriorExtension interior =
|
||
|
|
deformation::compileInteriorDeformationExtension(deformation::PowerLawRadialInteriorExtension{}, context);
|
||
|
|
const deformation::PreparedFixedInfinityRadialVacuumExtension vacuum =
|
||
|
|
deformation::compileVacuumDeformationExtension(deformation::FixedInfinityRadialVacuumExtension{}, context);
|
||
|
|
|
||
|
|
mfem::Vector surface(interior.surfaceDisplacementSize());
|
||
|
|
mfem::Vector direction(interior.surfaceDisplacementSize());
|
||
|
|
for (int dof = 0; dof < surface.Size(); ++dof) {
|
||
|
|
const double index = static_cast<double>(dof + 1);
|
||
|
|
surface(dof) = 0.01 * std::sin(0.17 * index);
|
||
|
|
direction(dof) = std::cos(0.13 * index) - 0.2 * std::sin(0.31 * index);
|
||
|
|
}
|
||
|
|
|
||
|
|
constexpr double step = 1.0e-6;
|
||
|
|
mfem::Vector plusSurface(surface);
|
||
|
|
mfem::Vector minusSurface(surface);
|
||
|
|
plusSurface.Add(step, direction);
|
||
|
|
minusSurface.Add(-step, direction);
|
||
|
|
|
||
|
|
auto checkLinearization = [&](const auto &prepared, const int volumeSize, const auto &build) {
|
||
|
|
mfem::Vector plus(volumeSize);
|
||
|
|
mfem::Vector minus(volumeSize);
|
||
|
|
mfem::Vector jacobian(volumeSize);
|
||
|
|
build(prepared, plusSurface, plus);
|
||
|
|
build(prepared, minusSurface, minus);
|
||
|
|
prepared.applyJacobian(surface, direction, jacobian);
|
||
|
|
|
||
|
|
mfem::Vector centeredDifference(plus);
|
||
|
|
centeredDifference -= minus;
|
||
|
|
centeredDifference /= 2.0 * step;
|
||
|
|
CHECK(radial_extension_test_utils::relativeError(jacobian, centeredDifference) < 2.0e-10);
|
||
|
|
|
||
|
|
mfem::Vector volumeDual(volumeSize);
|
||
|
|
for (int dof = 0; dof < volumeDual.Size(); ++dof) {
|
||
|
|
const double index = static_cast<double>(dof + 1);
|
||
|
|
volumeDual(dof) = std::sin(0.07 * index) + 0.3 * std::cos(0.11 * index);
|
||
|
|
}
|
||
|
|
mfem::Vector surfaceDual(surface.Size());
|
||
|
|
prepared.applyJacobianTranspose(surface, volumeDual, surfaceDual);
|
||
|
|
|
||
|
|
const double volumeWork =
|
||
|
|
radial_extension_test_utils::globalInnerProduct(jacobian, volumeDual, fem.mesh->GetComm());
|
||
|
|
const double surfaceWork =
|
||
|
|
radial_extension_test_utils::globalInnerProduct(direction, surfaceDual, fem.mesh->GetComm());
|
||
|
|
const double scale = std::max({1.0, std::abs(volumeWork), std::abs(surfaceWork)});
|
||
|
|
CHECK(std::abs(volumeWork - surfaceWork) <= 5.0e-13 * scale);
|
||
|
|
|
||
|
|
mfem::Vector pullback(surface.Size());
|
||
|
|
pullback = 1.0;
|
||
|
|
prepared.applyPullbackDerivative(surface, direction, volumeDual, pullback);
|
||
|
|
CHECK(pullback.Norml2() == 0.0);
|
||
|
|
};
|
||
|
|
|
||
|
|
checkLinearization(
|
||
|
|
interior, interior.interiorDisplacementSize(),
|
||
|
|
[](const auto &prepared, const mfem::Vector &input, mfem::Vector &output) {
|
||
|
|
prepared.buildInteriorDisplacement(input, output);
|
||
|
|
}
|
||
|
|
);
|
||
|
|
checkLinearization(
|
||
|
|
vacuum, vacuum.vacuumDisplacementSize(),
|
||
|
|
[](const auto &prepared, const mfem::Vector &input, mfem::Vector &output) {
|
||
|
|
prepared.buildVacuumDisplacement(input, output);
|
||
|
|
}
|
||
|
|
);
|
||
|
|
|
||
|
|
mfem::Vector wrongSurface(surface.Size() + 1);
|
||
|
|
mfem::Vector interiorOutput(interior.interiorDisplacementSize());
|
||
|
|
mfem::Vector vacuumOutput(vacuum.vacuumDisplacementSize());
|
||
|
|
CHECK_THROWS_AS(interior.buildInteriorDisplacement(wrongSurface, interiorOutput), std::invalid_argument);
|
||
|
|
CHECK_THROWS_AS(vacuum.buildVacuumDisplacement(wrongSurface, vacuumOutput), std::invalid_argument);
|
||
|
|
}
|
||
|
|
|
||
|
|
TEST_CASE(
|
||
|
|
"Logical Radial Deformation Remains Conforming And Orientation Preserving After Mesh Refinement",
|
||
|
|
tags::radial_deformation_extension_mapping
|
||
|
|
) {
|
||
|
|
namespace deformation = mean_field::deformation;
|
||
|
|
|
||
|
|
mean_field::utils::Args args = test_utils::setup_args();
|
||
|
|
mean_field::fem::FEM fem = mean_field::fem::setup_fem(args.mesh_file, args, 1);
|
||
|
|
REQUIRE(fem.okay());
|
||
|
|
|
||
|
|
deformation::PreparedNodalRadialSurface surface = radial_extension_test_utils::makePreparedSurface(fem);
|
||
|
|
const deformation::RadialDeformationExtensionCompilationContext context =
|
||
|
|
deformation::makeRadialDeformationExtensionCompilationContext<radial_extension_test_utils::Schema>(
|
||
|
|
*fem.surfaceDeformationFes, *fem.displacementFes, *fem.logicalReferenceMesh
|
||
|
|
);
|
||
|
|
deformation::PreparedPowerLawRadialInteriorExtension interior =
|
||
|
|
deformation::compileInteriorDeformationExtension(deformation::PowerLawRadialInteriorExtension{}, context);
|
||
|
|
deformation::PreparedFixedInfinityRadialVacuumExtension vacuum =
|
||
|
|
deformation::compileVacuumDeformationExtension(deformation::FixedInfinityRadialVacuumExtension{}, context);
|
||
|
|
auto prepared = deformation::composePreparedDomainDeformation(
|
||
|
|
std::move(surface), std::move(interior), std::move(vacuum), *fem.surfaceDeformationFes, *fem.displacementFes,
|
||
|
|
*fem.logicalReferenceMesh
|
||
|
|
);
|
||
|
|
|
||
|
|
mfem::Vector parameters(prepared.parameterCount());
|
||
|
|
mfem::Vector direction(prepared.parameterCount());
|
||
|
|
for (int parameter = 0; parameter < parameters.Size(); ++parameter) {
|
||
|
|
const double index = static_cast<double>(parameter + 1);
|
||
|
|
const double polarDirection = prepared.surfaceDeformationPrescription().radialDirection(parameter, 2);
|
||
|
|
const double quadrupoleValue = 0.5 * (3.0 * polarDirection * polarDirection - 1.0);
|
||
|
|
parameters(parameter) = 0.006 - 0.001 * quadrupoleValue;
|
||
|
|
direction(parameter) = std::sin(0.07 * index) - 0.4 * std::cos(0.13 * index);
|
||
|
|
}
|
||
|
|
|
||
|
|
mfem::Vector volumeDisplacement(prepared.volumeDisplacementSize());
|
||
|
|
prepared.buildVolumeDisplacement(parameters, volumeDisplacement);
|
||
|
|
const deformation::DomainDeformationGeometryReport geometry = prepared.inspectMappedGeometry(volumeDisplacement);
|
||
|
|
CAPTURE(geometry.minimumJacobianDeterminant);
|
||
|
|
REQUIRE(geometry.isOrientationPreserving());
|
||
|
|
|
||
|
|
mfem::Vector jacobianAction(prepared.volumeDisplacementSize());
|
||
|
|
mfem::Vector volumeDual(prepared.volumeDisplacementSize());
|
||
|
|
for (int dof = 0; dof < volumeDual.Size(); ++dof) {
|
||
|
|
const double index = static_cast<double>(dof + 1);
|
||
|
|
volumeDual(dof) = std::cos(0.017 * index) + 0.2 * std::sin(0.023 * index);
|
||
|
|
}
|
||
|
|
prepared.applyJacobian(parameters, direction, jacobianAction);
|
||
|
|
mfem::Vector parameterDual(prepared.parameterCount());
|
||
|
|
prepared.applyJacobianTranspose(parameters, volumeDual, parameterDual);
|
||
|
|
const double volumeWork =
|
||
|
|
radial_extension_test_utils::globalInnerProduct(jacobianAction, volumeDual, fem.mesh->GetComm());
|
||
|
|
const double parameterWork =
|
||
|
|
radial_extension_test_utils::globalInnerProduct(direction, parameterDual, fem.mesh->GetComm());
|
||
|
|
CHECK(
|
||
|
|
std::abs(volumeWork - parameterWork) <= 2.0e-12 * std::max({1.0, std::abs(volumeWork), std::abs(parameterWork)})
|
||
|
|
);
|
||
|
|
}
|