phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #4231 from mfem/2x2_matrix_exponential

Browse Source 
ExponentialMatrixCoefficient class
This commit is contained in:
Tzanio Kolev
2024-05-28 12:15:31 -07:00
committed by GitHub
parent 80f0f6cdb9 2636fffda9
commit 41a40ebf57
7 changed files with 395 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/bilininteg.hpp
+3 -3
View File
@@ -1741,7 +1741,7 @@ public:
{ vector_fe.CalcPhysDShape(Trans, shape); }
};
/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot \nabla \cdot u, \nabla \cdot v)$ in 2D
/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot \nabla u, \nabla \cdot v)$ in 2D
or 3D and where $\hat{V}$ is a vector coefficient, $u$ is in $H^1$ and $v$ is in $H(div)$. */
class MixedGradDivIntegrator : public MixedScalarVectorIntegrator
{
@@ -1780,7 +1780,7 @@ public:
{ scalar_fe.CalcPhysDivShape(Trans, shape); }
};
/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla \cdot v)$ in 2D
/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla v)$ in 2D
or 3D and where $\hat{V}$ is a vector coefficient, $u$ is in $H(div)$ and $v$ is in $H^1$. */
class MixedDivGradIntegrator : public MixedScalarVectorIntegrator
{
@@ -1820,7 +1820,7 @@ public:
{ scalar_fe.CalcPhysDivShape(Trans, shape); }
};
/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla \cdot v)$ in 2D or 3D
/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla v)$ in 2D or 3D
and where $\hat{V}$ is a vector coefficient, $u$ is in $H^1$ or $L_2$ and $v$ is in $H^1$. */
class MixedScalarWeakDivergenceIntegrator : public MixedScalarVectorIntegrator
{
fem/coefficient.cpp
+114
View File
@@ -924,6 +924,75 @@ void MatrixArrayCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
}
}
MatrixArrayVectorCoefficient::MatrixArrayVectorCoefficient (int dim)
: MatrixCoefficient (dim)
{
Coeff.SetSize(height);
ownCoeff.SetSize(height);
for (int i = 0; i < height; i++)
{
Coeff[i] = NULL;
ownCoeff[i] = true;
}
}
void MatrixArrayVectorCoefficient::SetTime(real_t t)
{
for (int i=0; i < height; i++)
{
if (Coeff[i]) { Coeff[i]->SetTime(t); }
}
this->MatrixCoefficient::SetTime(t);
}
void MatrixArrayVectorCoefficient::Set(int i, VectorCoefficient * c, bool own)
{
MFEM_ASSERT(i < height && i >= 0, "Row "
<< i << " does not exist. " <<
"Matrix height = " << height << ".");
if (ownCoeff[i]) { delete Coeff[i]; }
Coeff[i] = c;
ownCoeff[i] = own;
}
MatrixArrayVectorCoefficient::~MatrixArrayVectorCoefficient ()
{
for (int i=0; i < height; i++)
{
if (ownCoeff[i]) { delete Coeff[i]; }
}
}
void MatrixArrayVectorCoefficient::Eval(int i, Vector &V,
ElementTransformation &T,
const IntegrationPoint &ip)
{
MFEM_ASSERT(i < height && i >= 0, "Row "
<< i << " does not exist. " <<
"Matrix height = " << height << ".");
if (Coeff[i])
{
Coeff[i] -> Eval(V, T, ip);
}
else
{
V = 0.0;
}
}
void MatrixArrayVectorCoefficient::Eval(DenseMatrix &K,
ElementTransformation &T,
const IntegrationPoint &ip)
{
K.SetSize(height, width);
Vector V(width);
for (int i = 0; i < height; i++)
{
this->Eval(i, V, T, ip);
K.SetRow(i, V);
}
}
void MatrixRestrictedCoefficient::SetTime(real_t t)
{
if (c) { c->SetTime(t); }
@@ -1041,6 +1110,27 @@ real_t DeterminantCoefficient::Eval(ElementTransformation &T,
return ma.Det();
}
TraceCoefficient::TraceCoefficient(MatrixCoefficient &A)
: a(&A), ma(A.GetHeight(), A.GetWidth())
{
MFEM_ASSERT(A.GetHeight() == A.GetWidth(),
"TraceCoefficient: "
"Argument must be a square matrix.");
}
void TraceCoefficient::SetTime(real_t t)
{
if (a) { a->SetTime(t); }
this->Coefficient::SetTime(t);
}
real_t TraceCoefficient::Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{
a->Eval(ma, T, ip);
return ma.Trace();
}
VectorSumCoefficient::VectorSumCoefficient(int dim)
: VectorCoefficient(dim),
ACoef(NULL), BCoef(NULL),
@@ -1326,6 +1416,30 @@ void InverseMatrixCoefficient::Eval(DenseMatrix &M,
M.Invert();
}
ExponentialMatrixCoefficient::ExponentialMatrixCoefficient(MatrixCoefficient &A)
: MatrixCoefficient(A.GetHeight(), A.GetWidth()), a(&A)
{
MFEM_ASSERT(A.GetHeight() == A.GetWidth() && A.GetHeight() == 2,
"ExponentialMatrixCoefficient: "
<< "Argument must be a square 2x2 matrix."
<< " Height = " << A.GetHeight()
<< ", Width = " << A.GetWidth());
}
void ExponentialMatrixCoefficient::SetTime(real_t t)
{
if (a) { a->SetTime(t); }
this->MatrixCoefficient::SetTime(t);
}
void ExponentialMatrixCoefficient::Eval(DenseMatrix &M,
ElementTransformation &T,
const IntegrationPoint &ip)
{
a->Eval(M, T, ip);
M.Exponential();
}
OuterProductCoefficient::OuterProductCoefficient(VectorCoefficient &A,
VectorCoefficient &B)
: MatrixCoefficient(A.GetVDim(), B.GetVDim()), a(&A), b(&B),
fem/coefficient.hpp
+90 -2
View File
@@ -1334,6 +1334,46 @@ public:
virtual ~MatrixArrayCoefficient();
};
/** @brief Matrix coefficient defined row-wise by an array of vector
coefficients. Rows that are not set will evaluate to zero. The
matrix coefficient is stored as an array indexing the rows of
the matrix. */
class MatrixArrayVectorCoefficient : public MatrixCoefficient
{
private:
Array<VectorCoefficient *> Coeff;
Array<bool> ownCoeff;
public:
/** @brief Construct a coefficient matrix of dimensions @a dim * @a dim. The
actual coefficients still need to be added with Set(). */
explicit MatrixArrayVectorCoefficient (int dim);
/// Set the time for internally stored coefficients
void SetTime(real_t t) override;
/// Get the vector coefficient located at the i-th row of the matrix
VectorCoefficient* GetCoeff (int i) { return Coeff[i]; }
/** @brief Set the coefficient located at the i-th row of the matrix.
By this will take ownership of the Coefficient passed in, but this
can be overridden with the @a own parameter. */
void Set(int i, VectorCoefficient * c, bool own=true);
using MatrixCoefficient::Eval;
/// Evaluate coefficient located at the i-th row of the matrix using integration
/// point @a ip.
void Eval(int i, Vector &V, ElementTransformation &T,
const IntegrationPoint &ip);
/// Evaluate the matrix coefficient @a ip.
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip) override;
virtual ~MatrixArrayVectorCoefficient();
};
/** @brief Derived matrix coefficient that has the value of the parent matrix
coefficient where it is active and is zero otherwise. */
@@ -1767,6 +1807,31 @@ public:
const IntegrationPoint &ip);
};
/// Scalar coefficient defined as the trace of a matrix coefficient
class TraceCoefficient : public Coefficient
{
private:
MatrixCoefficient * a;
mutable DenseMatrix ma;
public:
/// Construct with the matrix.
TraceCoefficient(MatrixCoefficient &A);
/// Set the time for internally stored coefficients
void SetTime(real_t t);
/// Reset the matrix coefficient
void SetACoef(MatrixCoefficient &A) { a = &A; }
/// Return the matrix coefficient
MatrixCoefficient * GetACoef() const { return a; }
/// Evaluate the trace coefficient at @a ip.
virtual real_t Eval(ElementTransformation &T,
const IntegrationPoint &ip);
};
/// Vector coefficient defined as the linear combination of two vectors
class VectorSumCoefficient : public VectorCoefficient
{
@@ -2118,7 +2183,7 @@ public:
const IntegrationPoint &ip);
};
/// Matrix coefficient defined as the transpose a matrix coefficient
/// Matrix coefficient defined as the transpose of a matrix coefficient
class TransposeMatrixCoefficient : public MatrixCoefficient
{
private:
@@ -2141,7 +2206,7 @@ public:
const IntegrationPoint &ip);
};
/// Matrix coefficient defined as the inverse a matrix coefficient.
/// Matrix coefficient defined as the inverse of a matrix coefficient.
class InverseMatrixCoefficient : public MatrixCoefficient
{
private:
@@ -2164,6 +2229,29 @@ public:
const IntegrationPoint &ip);
};
/// Matrix coefficient defined as the exponential of a matrix coefficient.
class ExponentialMatrixCoefficient : public MatrixCoefficient
{
private:
MatrixCoefficient * a;
public:
/// Construct the matrix coefficient. Result is $ \exp(A) $.
ExponentialMatrixCoefficient(MatrixCoefficient &A);
/// Set the time for internally stored coefficients
void SetTime(real_t t);
/// Reset the matrix coefficient
void SetACoef(MatrixCoefficient &A) { a = &A; }
/// Return the matrix coefficient
MatrixCoefficient * GetACoef() const { return a; }
/// Evaluate the matrix coefficient at @a ip.
virtual void Eval(DenseMatrix &M, ElementTransformation &T,
const IntegrationPoint &ip);
};
/// Matrix coefficient defined as the outer product of two vector coefficients.
class OuterProductCoefficient : public MatrixCoefficient
{
linalg/densemat.cpp
+63
View File
@@ -532,6 +532,69 @@ MatrixInverse *DenseMatrix::Inverse() const
return new DenseMatrixInverse(*this);
}
void DenseMatrix::Exponential()
{
MFEM_ASSERT(Height() == Width() && Height() <= 2,
"The matrix must be square and "
<< "of size less than or equal to 2."
<< " Height() = " << Height()
<< ", Width() = " << Width());
switch (Height())
{
case 1:
{
data[0] = std::exp(data[0]);
break;
}
case 2:
{
/// Formulas from Corollary 2.4 of doi:10.1109/9.233156
/// Note typo in the paper, in the prefactor in the equation under (i).
const real_t a = data[0];
const real_t b = data[1];
const real_t c = data[2];
const real_t d = data[3];
const real_t e = (a - d)*(a - d) + 4*b*c;
const real_t f = std::exp((a + d)/2.0);
const real_t g = std::sqrt(std::abs(e)) / 2.0;
if (e == 0)
{
data[0] = 1.0 + (a - d)/2.0;
data[3] = 1.0 - (a - d)/2.0;
}
else if (e > 0)
{
data[0] = std::cosh(g) + (a - d)/2 * std::sinh(g) / g;
data[1] = b * std::sinh(g) / g;
data[2] = c * std::sinh(g) / g;
data[3] = std::cosh(g) - (a - d)/2 * std::sinh(g) / g;
}
else
{
data[0] = std::cos(g) + (a - d)/2 * std::sin(g) / g;
data[1] = b * std::sin(g) / g;
data[2] = c * std::sin(g) / g;
data[3] = std::cos(g) - (a - d)/2 * std::sin(g) / g;
}
for (int i = 0; i < 4; i++)
{
data[i] *= f;
}
break;
}
case 3:
{
MFEM_ABORT("3x3 matrices are not currently supported");
}
default:
{
MFEM_ABORT("Only 1x1 and 2x2 matrices are currently supported");
}
}
}
real_t DenseMatrix::Det() const
{
MFEM_ASSERT(Height() == Width() && Height() > 0,
linalg/densemat.hpp
+4
View File
@@ -207,6 +207,10 @@ public:
/// Replaces the current matrix with its square root inverse
void SquareRootInverse();
/// Replaces the current matrix with its exponential
/// (currently only supports 2x2 matrices)
void Exponential();
/// Calculates the determinant of the matrix
/// (optimized for 2x2, 3x3, and 4x4 matrices)
real_t Det() const;
tests/unit/fem/test_coefficient.cpp
+60
View File
@@ -304,6 +304,66 @@ TEST_CASE("Piecewise Matrix Coefficient", "[Coefficient]")
}
}
TEST_CASE("MatrixArrayVectorCoefficient", "[Coefficient]")
{
Vector V1(2), V2(2);
V1(0) = 0.0; V1(1) = 1.0;
V2(0) = 2.0; V2(1) = 3.0;
VectorConstantCoefficient Coef1(V1), Coef2(V2);
IsoparametricTransformation T;
IntegrationPoint ip;
MatrixArrayVectorCoefficient mavc(2);
Vector V(2);
// Verify zeros for unset rows
int row = 0;
mavc.Eval(row, V, T, ip);
REQUIRE(V(0) == MFEM_Approx(0.0));
REQUIRE(V(1) == MFEM_Approx(0.0));
row = 1;
mavc.Eval(row, V, T, ip);
REQUIRE(V(0) == MFEM_Approx(0.0));
REQUIRE(V(1) == MFEM_Approx(0.0));
DenseMatrix K(2);
mavc.Eval(K, T, ip);
REQUIRE(K(0,0) == MFEM_Approx(0.0));
REQUIRE(K(0,1) == MFEM_Approx(0.0));
REQUIRE(K(1,0) == MFEM_Approx(0.0));
REQUIRE(K(1,1) == MFEM_Approx(0.0));
// Test setting individual rows
row = 0;
mavc.Set(row, &Coef1, false);
mavc.Eval(row, V, T, ip);
REQUIRE(V(0) == MFEM_Approx(0.0));
REQUIRE(V(1) == MFEM_Approx(1.0));
row = 1;
mavc.Eval(row, V, T, ip);
REQUIRE(V(0) == MFEM_Approx(0.0));
REQUIRE(V(1) == MFEM_Approx(0.0));
mavc.Set(row, &Coef2, false);
row = 0;
mavc.Eval(row, V, T, ip);
REQUIRE(V(0) == MFEM_Approx(0.0));
REQUIRE(V(1) == MFEM_Approx(1.0));
row = 1;
mavc.Eval(row, V, T, ip);
REQUIRE(V(0) == MFEM_Approx(2.0));
REQUIRE(V(1) == MFEM_Approx(3.0));
mavc.Eval(K, T, ip);
REQUIRE(K(0,0) == MFEM_Approx(0.0));
REQUIRE(K(0,1) == MFEM_Approx(1.0));
REQUIRE(K(1,0) == MFEM_Approx(2.0));
REQUIRE(K(1,1) == MFEM_Approx(3.0));
}
TEST_CASE("Symmetric Matrix Coefficient", "[Coefficient]")
{
int d = 3;
tests/unit/linalg/test_matrix_dense.cpp
+61
View File
@@ -623,6 +623,67 @@ TEST_CASE("MatrixInverse", "[DenseMatrix]")
}
}
TEST_CASE("Exponential", "[DenseMatrix]")
{
// case 1
DenseMatrix A(2,2);
A(0,0) = 5.0;
A(0,1) = 3.0;
A(1,0) = 0.0;
A(1,1) = 5.0;
A.Exponential();
DenseMatrix expA(2,2);
expA(0,0) = std::exp(5.0);
expA(0,1) = 3.0 * std::exp(5.0);
expA(1,0) = 0.0;
expA(1,1) = std::exp(5.0);
A.Print();
expA.Print();
REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
// case 2
A(0,0) = 3.0;
A(0,1) = 5.0;
A(1,0) = 4.0;
A(1,1) = 2.0;
A.Exponential();
expA(0,0) = 4.0 / (9.0 * std::exp(2.0)) + (5.0 * std::exp(7.0)) / 9.0;
expA(0,1) = (5.0 * std::exp(7.0)) / 9.0 - 5.0 / (9.0 * std::exp(2.0));
expA(1,0) = (4.0 * std::exp(7.0)) / 9.0 - 4.0 / (9.0 * std::exp(2.0));
expA(1,1) = 5.0 / (9.0 * std::exp(2.0)) + (4.0 * std::exp(7.0)) / 9.0;
REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
// case 3
A(0,0) = 10.0;
A(0,1) = 2.0;
A(1,0) = -2.0;
A(1,1) = 8.0;
A.Exponential();
expA(0,0) = std::exp(9.0) * (std::sin(std::sqrt(3.0)) / std::sqrt(3.0)
+ std::cos(std::sqrt(3.0)));
expA(0,1) = 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
expA(1,0) = - 2.0 * std::exp(9.0) * std::sin(std::sqrt(3.0)) / std::sqrt(3.0);
expA(1,1) = std::exp(9.0) * (std::cos(std::sqrt(3.0))
- std::sin(std::sqrt(3.0)) / std::sqrt(3.0));
REQUIRE(A(0,0) == MFEM_Approx(expA(0,0)));
REQUIRE(A(0,1) == MFEM_Approx(expA(0,1)));
REQUIRE(A(1,0) == MFEM_Approx(expA(1,0)));
REQUIRE(A(1,1) == MFEM_Approx(expA(1,1)));
}
#ifdef MFEM_USE_LAPACK
enum class TestCase { GenEigSPD, GenEigGE, SVD};