feat(mfem-backend): ported DofMap system from MeanField
Broadly MeanField had a system that would map the full index of degrees of freedom to some reduced index space based on domain support. This system has been cleaned up and ported. Further, we have consolidated this system into the discretization system from where it lived in MeanField (in the operators subsystem)
This commit is contained in:
@@ -0,0 +1,174 @@
|
||||
#include "mfem/fem/fe/fe_base.hpp"
|
||||
#include "mfem/fem/intrules.hpp"
|
||||
#include "mfem/linalg/vector.hpp"
|
||||
#include "mfem/linalg/densemat.hpp"
|
||||
|
||||
#include "serif/discretization/backend/mfem/basis_evaluation.hpp"
|
||||
|
||||
#include "serif/discretization/backend/mfem/tensor_basis.hpp"
|
||||
#include "serif/utils/error/errors.hpp"
|
||||
|
||||
namespace serif::discretization::backend::mfem {
|
||||
void BasisEvaluation::sample(
|
||||
const ::mfem::FiniteElement& element,
|
||||
const ::mfem::IntegrationRule& rule,
|
||||
const ReferenceOperation operation,
|
||||
const int copies,
|
||||
TensorBasisCache& cache,
|
||||
const ContractionMode mode
|
||||
) {
|
||||
using namespace utils::errors;
|
||||
|
||||
const bool vector = element.GetRangeType() == ::mfem::FiniteElement::VECTOR;
|
||||
const bool operation_supported = vector ? (operation == ReferenceOperation::value || operation == ReferenceOperation::divergence) : operation != ReferenceOperation::divergence;
|
||||
const bool mapping_supported = vector ? copies == 1 && element.GetMapType() == ::mfem::FiniteElement::H_DIV : element.GetMapType() == ::mfem::FiniteElement::VALUE;
|
||||
|
||||
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(element.GetDim() == 3 && copies > 0 && operation_supported && mapping_supported, "Canonical basis sampling requires 3D scalar VALUE bases or intrinsic H(div) bases and a compatible reference operation.");
|
||||
|
||||
const int samples = operation == ReferenceOperation::gradient ? 3 : operation == ReferenceOperation::value_and_gradient ? 4 : 1;
|
||||
m_components = vector ? (operation == ReferenceOperation::value ? 3 : 1) : copies * samples;
|
||||
m_points = rule.GetNPoints();
|
||||
m_width = copies * element.GetDof();
|
||||
m_tensor = mode == ContractionMode::automatic ? cache.get(element, rule, operation, copies) : nullptr;
|
||||
|
||||
if (m_tensor) {
|
||||
m_dense.SetSize(0, 0);
|
||||
m_mapping.SetSize(m_points * m_components, m_components);
|
||||
m_mapping = 0.0;
|
||||
|
||||
for (int q = 0; q < m_points; ++q) {
|
||||
for (int component = 0; component < m_components; ++component) {
|
||||
m_mapping(q * m_components + component, component) = 1.0;
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
m_mapping.SetSize(0, 0);
|
||||
m_dense.SetSize(m_points * m_components, m_width);
|
||||
m_dense = 0.0;
|
||||
|
||||
::mfem::Vector shape(element.GetDof());
|
||||
::mfem::DenseMatrix gradients(element.GetDof(), 3);
|
||||
|
||||
for (int q = 0; q < m_points; ++q) {
|
||||
if (vector) {
|
||||
if (operation == ReferenceOperation::divergence) {
|
||||
element.CalcDivShape(rule.IntPoint(q), shape);
|
||||
} else {
|
||||
element.CalcVShape(rule.IntPoint(q), gradients);
|
||||
}
|
||||
} else {
|
||||
element.CalcShape(rule.IntPoint(q), shape);
|
||||
if (samples > 1) {
|
||||
element.CalcDShape(rule.IntPoint(q), gradients);
|
||||
}
|
||||
}
|
||||
|
||||
for (int copy = 0; copy < copies; ++copy) {
|
||||
for (int sample = 0; sample < (vector ? m_components : samples); ++sample) {
|
||||
const int derivative = operation == ReferenceOperation::gradient ? sample : sample - 1;
|
||||
|
||||
for (int dof = 0; dof < element.GetDof(); ++dof) {
|
||||
double value;
|
||||
|
||||
if (vector) {
|
||||
value = operation == ReferenceOperation::divergence ? shape[dof] : gradients(dof, sample);
|
||||
} else {
|
||||
value = samples == 1 || derivative < 0 ? shape[dof] : gradients(dof, derivative);
|
||||
}
|
||||
|
||||
m_dense(q * m_components + copy * samples + sample, copy * element.GetDof() + dof) = value;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void BasisEvaluation::set_constant(const int points) {
|
||||
m_tensor.reset();
|
||||
m_components = 1;
|
||||
m_points = points;
|
||||
m_width = 1;
|
||||
m_dense.SetSize(points, 1);
|
||||
m_dense = 1.0;
|
||||
m_mapping.SetSize(0, 0);
|
||||
}
|
||||
|
||||
|
||||
void BasisEvaluation::reset_from(const BasisEvaluation& reference, const int components) {
|
||||
m_tensor = reference.m_tensor;
|
||||
m_points = reference.m_points;
|
||||
m_width = reference.m_width;
|
||||
m_components = components;
|
||||
|
||||
if (m_tensor) {
|
||||
m_mapping.SetSize(m_points * components, m_tensor->components());
|
||||
m_dense.SetSize(0, 0);
|
||||
} else {
|
||||
m_dense.SetSize(m_points * components, m_width);
|
||||
m_mapping.SetSize(0, 0);
|
||||
}
|
||||
}
|
||||
|
||||
void BasisEvaluation::evaluate(const ::mfem::Vector& coefficients, ::mfem::Vector& values, BasisWorkspace& workspace) const {
|
||||
using namespace utils::errors;
|
||||
|
||||
serif_ensure<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(
|
||||
coefficients.Size() == m_width,
|
||||
"Basis evaluation coefficient count does not match its element."
|
||||
);
|
||||
values.SetSize(height());
|
||||
if (!m_tensor) {
|
||||
m_dense.Mult(coefficients, values);
|
||||
return;
|
||||
}
|
||||
|
||||
m_tensor->apply(coefficients, workspace.reference_values, workspace.tensor);
|
||||
values = 0.0;
|
||||
for (int q = 0; q < m_points; ++q) {
|
||||
for (int output = 0; output < m_components; ++output) {
|
||||
for (int input = 0; input < m_tensor->components(); ++input) {
|
||||
values[q * m_components + output] += m_mapping(q * m_components + output, input) * workspace.reference_values[q * m_tensor->components() + input];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void BasisEvaluation::integrate(const ::mfem::Vector& values, ::mfem::Vector& coefficients, BasisWorkspace& workspace) const {
|
||||
using namespace utils::errors;
|
||||
|
||||
serif_ensure<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(
|
||||
values.Size() == height(),
|
||||
"Basis transpose quadrature size does not match its rule."
|
||||
);
|
||||
coefficients.SetSize(m_width);
|
||||
if (!m_tensor) {
|
||||
m_dense.MultTranspose(values, coefficients);
|
||||
return;
|
||||
}
|
||||
|
||||
workspace.reference_values.SetSize(m_points * m_tensor->components());
|
||||
workspace.reference_values = 0.0;
|
||||
for (int q = 0; q < m_points; ++q) {
|
||||
for (int input = 0; input < m_tensor->components(); ++input) {
|
||||
for (int output = 0; output < m_components; ++output) {
|
||||
workspace.reference_values[q * m_tensor->components() + input] += m_mapping(q * m_components + output, input) * values[q * m_components + output];
|
||||
}
|
||||
}
|
||||
}
|
||||
m_tensor->apply(workspace.reference_values, coefficients, workspace.tensor, true);
|
||||
}
|
||||
|
||||
int BasisEvaluation::height() const noexcept {
|
||||
return m_points * m_components;
|
||||
}
|
||||
|
||||
int BasisEvaluation::width() const noexcept {
|
||||
return m_width;
|
||||
}
|
||||
|
||||
bool BasisEvaluation::uses_tensor() const noexcept {
|
||||
return static_cast<bool>(m_tensor);
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user