15#ifndef PBAT_FEM_MASS_H
16#define PBAT_FEM_MASS_H
26#include <tbb/parallel_for.h>
41template <
typename TElement,
typename TN>
43ElementMassMatrix(Eigen::MatrixBase<TN>
const& N,
typename TN::Scalar w,
typename TN::Scalar rho)
46 TElement::kNodes == TN::RowsAtCompileTime || TN::RowsAtCompileTime == Eigen::Dynamic,
47 "Shape function vector size must match number of element nodes");
48 return w * rho * N * N.transpose();
69 typename TDerivedrhog,
72 Eigen::MatrixBase<TN>
const& Neg,
73 Eigen::DenseBase<TDerivedwg>
const& wg,
74 Eigen::DenseBase<TDerivedrhog>
const& rhog,
75 Eigen::PlainObjectBase<TDerivedOut>& Meg)
77 constexpr int kNodes = TElement::kNodes;
78 const auto nQuadPts = Neg.cols();
79 Meg.resize(kNodes, kNodes * nQuadPts);
80 for (Eigen::Index g = 0; g < nQuadPts; ++g)
82 auto const N = Neg.col(g);
84 auto const rho = rhog(g);
85 auto Me = w * rho * N * N.transpose();
86 Meg.template block<kNodes, kNodes>(0, g * kNodes) = Me;
102template <
typename TElement,
typename TN,
typename TDerivedwg,
typename TDerivedrhog>
104 Eigen::MatrixBase<TN>
const& Neg,
105 Eigen::DenseBase<TDerivedwg>
const& wg,
106 Eigen::DenseBase<TDerivedrhog>
const& rhog)
108 using ScalarType =
typename TN::Scalar;
109 constexpr int kNodes = TElement::kNodes;
110 Eigen::Matrix<ScalarType, kNodes, Eigen::Dynamic> Meg;
144 Eigen::DenseBase<TDerivedE>
const& E,
146 Eigen::DenseBase<TDerivedeg>
const& eg,
147 Eigen::MatrixBase<TDerivedMe>
const& Me,
149 Eigen::MatrixBase<TDerivedIn>
const& X,
150 Eigen::DenseBase<TDerivedOut>& Y);
168template <CMesh TMesh,
class TDerivedeg,
class TDerivedMe,
class TDerivedIn,
class TDerivedOut>
171 Eigen::DenseBase<TDerivedeg>
const& eg,
172 Eigen::MatrixBase<TDerivedMe>
const& Me,
174 Eigen::MatrixBase<TDerivedIn>
const& X,
175 Eigen::DenseBase<TDerivedOut>& Y)
179 static_cast<typename TMesh::IndexType
>(mesh.X.cols()),
209 Eigen::StorageOptions Options,
216 Eigen::DenseBase<TDerivedE>
const& E,
218 Eigen::DenseBase<TDerivedeg>
const& eg,
219 Eigen::DenseBase<TDerivedwg>
const& wg,
220 Eigen::DenseBase<TDerivedrhog>
const& rhog,
221 Eigen::MatrixBase<TDerivedNeg>
const& Neg,
223 -> Eigen::SparseMatrix<typename TDerivedNeg::Scalar, Options, typename TDerivedE::Scalar>;
245 Eigen::StorageOptions Options,
250 Eigen::DenseBase<TDerivedE>
const& E,
252 Eigen::DenseBase<TDerivedeg>
const& eg,
253 Eigen::MatrixBase<TDerivedMe>
const& Meg,
255 -> Eigen::SparseMatrix<typename TDerivedMe::Scalar, Options, typename TDerivedE::Scalar>;
274 Eigen::StorageOptions Options,
282 Eigen::DenseBase<TDerivedeg>
const& eg,
283 Eigen::DenseBase<TDerivedwg>
const& wg,
284 Eigen::DenseBase<TDerivedrhog>
const& rhog,
285 Eigen::MatrixBase<TDerivedNeg>
const& Neg,
287 -> Eigen::SparseMatrix<typename TDerivedNeg::Scalar, Options, typename TMesh::IndexType>
291 static_cast<typename TMesh::IndexType
>(mesh.X.cols()),
313template <Eigen::StorageOptions Options, CMesh TMesh,
class TDerivedeg,
class TDerivedMe>
316 Eigen::DenseBase<TDerivedeg>
const& eg,
317 Eigen::MatrixBase<TDerivedMe>
const& Meg,
319 -> Eigen::SparseMatrix<typename TDerivedMe::Scalar, Options, typename TMesh::IndexType>
323 static_cast<typename TMesh::IndexType
>(mesh.X.cols()),
358 Eigen::DenseBase<TDerivedE>
const& E,
360 Eigen::DenseBase<TDerivedeg>
const& eg,
361 Eigen::DenseBase<TDerivedwg>
const& wg,
362 Eigen::DenseBase<TDerivedrhog>
const& rhog,
363 Eigen::MatrixBase<TDerivedNeg>
const& Neg,
365 Eigen::PlainObjectBase<TDerivedOut>& m);
392 Eigen::DenseBase<TDerivedE>
const& E,
394 Eigen::DenseBase<TDerivedeg>
const& eg,
395 Eigen::MatrixBase<TDerivedMe>
const& Meg,
397 Eigen::PlainObjectBase<TDerivedOut>& m);
425 Eigen::DenseBase<TDerivedeg>
const& eg,
426 Eigen::DenseBase<TDerivedwg>
const& wg,
427 Eigen::DenseBase<TDerivedrhog>
const& rhog,
428 Eigen::MatrixBase<TDerivedNeg>
const& Neg,
430 Eigen::PlainObjectBase<TDerivedOut>& m)
434 static_cast<typename TMesh::IndexType
>(mesh.X.cols()),
457template <CMesh TMesh,
class TDerivedeg,
class TDerivedMe,
class TDerivedOut>
460 Eigen::DenseBase<TDerivedeg>
const& eg,
461 Eigen::MatrixBase<TDerivedMe>
const& Meg,
463 Eigen::PlainObjectBase<TDerivedOut>& m)
467 static_cast<typename TMesh::IndexType
>(mesh.X.cols()),
501 Eigen::DenseBase<TDerivedE>
const& E,
503 Eigen::DenseBase<TDerivedeg>
const& eg,
504 Eigen::DenseBase<TDerivedwg>
const& wg,
505 Eigen::DenseBase<TDerivedrhog>
const& rhog,
506 Eigen::MatrixBase<TDerivedNeg>
const& Neg,
507 int dims = 1) -> Eigen::Vector<typename TDerivedNeg::Scalar, Eigen::Dynamic>
509 using ScalarType =
typename TDerivedNeg::Scalar;
510 auto const numberOfDofs = dims * nNodes;
511 Eigen::Vector<ScalarType, Eigen::Dynamic> m(numberOfDofs);
532template <CElement TElement,
int Dims,
class TDerivedE,
class TDerivedeg,
class TDerivedMe>
534 Eigen::DenseBase<TDerivedE>
const& E,
536 Eigen::DenseBase<TDerivedeg>
const& eg,
537 Eigen::MatrixBase<TDerivedMe>
const& Meg,
538 int dims = 1) -> Eigen::Vector<typename TDerivedMe::Scalar, Eigen::Dynamic>
540 using ScalarType =
typename TDerivedMe::Scalar;
541 auto const numberOfDofs = dims * nNodes;
542 Eigen::Vector<ScalarType, Eigen::Dynamic> m(numberOfDofs);
550template <CMesh TMesh,
class TDerivedeg,
class TDerivedwg,
class TDerivedrhog,
class TDerivedNeg>
553 Eigen::DenseBase<TDerivedeg>
const& eg,
554 Eigen::DenseBase<TDerivedwg>
const& wg,
555 Eigen::DenseBase<TDerivedrhog>
const& rhog,
556 Eigen::MatrixBase<TDerivedNeg>
const& Neg,
561 static_cast<typename TMesh::IndexType
>(mesh.X.cols()),
573template <CMesh TMesh,
class TDerivedeg,
class TDerivedMe>
576 Eigen::DenseBase<TDerivedeg>
const& eg,
577 Eigen::MatrixBase<TDerivedMe>
const& Meg,
582 static_cast<typename TMesh::IndexType
>(mesh.X.cols()),
598 Eigen::DenseBase<TDerivedE>
const& E,
600 Eigen::DenseBase<TDerivedeg>
const& eg,
601 Eigen::MatrixBase<TDerivedMe>
const& Meg,
603 Eigen::MatrixBase<TDerivedIn>
const& X,
604 Eigen::DenseBase<TDerivedOut>& Y)
606 PBAT_PROFILE_NAMED_SCOPE(
"pbat.fem.GemmMass");
609 auto const numberOfDofs = dims * nNodes;
610 bool const bDimensionsMatch =
611 (X.cols() == Y.cols()) and (X.rows() == numberOfDofs) and (Y.rows() == numberOfDofs);
612 if (not bDimensionsMatch)
614 std::string
const what = fmt::format(
615 "Expected input and output to have {} rows and same number of columns, but got "
616 "dimensions X,Y=({} x {}), ({} x {})",
622 throw std::invalid_argument(what);
627 std::string
const what = fmt::format(
"Expected dims >= 1, got {} instead", dims);
628 throw std::invalid_argument(what);
631 auto constexpr kNodesPerElement = TElement::kNodes;
632 auto const nQuadPts = eg.size();
635 for (
auto c = 0; c < X.cols(); ++c)
637 for (
auto g = 0; g < nQuadPts; ++g)
639 auto const e = eg(g);
640 auto const nodes = E.col(e);
644 Meg.template block<kNodesPerElement, kNodesPerElement>(0, g * kNodesPerElement);
647 auto ye = Y.col(c).reshaped(dims, nNodes)(Eigen::placeholders::all, nodes);
648 auto const xe = X.col(c).reshaped(dims, nNodes)(Eigen::placeholders::all, nodes);
649 ye += xe * Me.transpose();
658 Eigen::StorageOptions Options,
665 Eigen::DenseBase<TDerivedE>
const& E,
667 Eigen::DenseBase<TDerivedeg>
const& eg,
668 Eigen::DenseBase<TDerivedwg>
const& wg,
669 Eigen::DenseBase<TDerivedrhog>
const& rhog,
670 Eigen::MatrixBase<TDerivedNeg>
const& Neg,
672 -> Eigen::SparseMatrix<typename TDerivedNeg::Scalar, Options, typename TDerivedE::Scalar>
674 PBAT_PROFILE_NAMED_SCOPE(
"pbat.fem.MassMatrix");
675 using ScalarType =
typename TDerivedNeg::Scalar;
676 using IndexType =
typename TDerivedE::Scalar;
677 using SparseMatrixType = Eigen::SparseMatrix<ScalarType, Options, IndexType>;
678 using Triplet = Eigen::Triplet<ScalarType, IndexType>;
680 auto constexpr kNodesPerElement = TElement::kNodes;
681 auto const nQuadPts = eg.size();
683 std::vector<Triplet> triplets{};
685 static_cast<std::size_t
>(kNodesPerElement * kNodesPerElement * nQuadPts * dims));
687 auto const numberOfDofs = dims * nNodes;
688 SparseMatrixType M(numberOfDofs, numberOfDofs);
690 for (
auto g = 0; g < nQuadPts; ++g)
692 auto const e = eg(g);
693 auto const nodes = E.col(e);
694 auto const w = wg(g);
695 auto const rho = rhog(g);
698 auto const N = Neg.col(g);
701 auto const Me = w * rho * N * N.transpose();
704 for (
auto i = 0; i < kNodesPerElement; ++i)
706 for (
auto j = 0; j < kNodesPerElement; ++j)
708 for (
auto d = 0; d < dims; ++d)
710 auto const ni =
static_cast<IndexType
>(dims * nodes(i) + d);
711 auto const nj =
static_cast<IndexType
>(dims * nodes(j) + d);
712 triplets.emplace_back(ni, nj, Me(i, j));
718 M.setFromTriplets(triplets.begin(), triplets.end());
725 Eigen::StorageOptions Options,
730 Eigen::DenseBase<TDerivedE>
const& E,
732 Eigen::DenseBase<TDerivedeg>
const& eg,
733 Eigen::MatrixBase<TDerivedMe>
const& Meg,
735 -> Eigen::SparseMatrix<typename TDerivedMe::Scalar, Options, typename TDerivedE::Scalar>
737 PBAT_PROFILE_NAMED_SCOPE(
"pbat.fem.MassMatrix");
738 using ScalarType =
typename TDerivedMe::Scalar;
739 using IndexType =
typename TDerivedE::Scalar;
740 using SparseMatrixType = Eigen::SparseMatrix<ScalarType, Options, IndexType>;
741 using Triplet = Eigen::Triplet<ScalarType, IndexType>;
743 auto constexpr kNodesPerElement = TElement::kNodes;
744 auto const nQuadPts = eg.size();
746 std::vector<Triplet> triplets{};
748 static_cast<std::size_t
>(kNodesPerElement * kNodesPerElement * nQuadPts * dims));
750 auto const numberOfDofs = dims * nNodes;
751 SparseMatrixType M(numberOfDofs, numberOfDofs);
753 for (
auto g = 0; g < nQuadPts; ++g)
755 auto const e = eg(g);
756 auto const nodes = E.col(e);
760 Meg.template block<kNodesPerElement, kNodesPerElement>(0, g * kNodesPerElement);
763 for (
auto i = 0; i < kNodesPerElement; ++i)
765 for (
auto j = 0; j < kNodesPerElement; ++j)
767 for (
auto d = 0; d < dims; ++d)
769 auto const ni =
static_cast<IndexType
>(dims * nodes(i) + d);
770 auto const nj =
static_cast<IndexType
>(dims * nodes(j) + d);
771 triplets.emplace_back(ni, nj, Me(i, j));
777 M.setFromTriplets(triplets.begin(), triplets.end());
792 Eigen::DenseBase<TDerivedE>
const& E,
794 Eigen::DenseBase<TDerivedeg>
const& eg,
795 Eigen::DenseBase<TDerivedwg>
const& wg,
796 Eigen::DenseBase<TDerivedrhog>
const& rhog,
797 Eigen::MatrixBase<TDerivedNeg>
const& Neg,
799 Eigen::PlainObjectBase<TDerivedOut>& m)
801 auto const numberOfDofs = dims * nNodes;
802 m.resize(numberOfDofs, 1);
804 auto constexpr kNodesPerElement = TElement::kNodes;
805 auto const nQuadPts = eg.size();
806 for (
auto g = 0; g < nQuadPts; ++g)
808 auto const e = eg(g);
809 auto const nodes = E.col(e);
810 auto const w = wg(g);
811 auto const rho = rhog(g);
814 auto const N = Neg.col(g);
817 auto const Me = (w * rho * N * N.transpose()).eval();
820 for (
auto i = 0; i < kNodesPerElement; ++i)
822 for (
auto j = 0; j < kNodesPerElement; ++j)
824 for (
auto d = 0; d < dims; ++d)
826 auto const ni = dims * nodes(i) + d;
842 Eigen::DenseBase<TDerivedE>
const& E,
844 Eigen::DenseBase<TDerivedeg>
const& eg,
845 Eigen::MatrixBase<TDerivedMe>
const& Meg,
847 Eigen::PlainObjectBase<TDerivedOut>& m)
849 auto const numberOfDofs = dims * nNodes;
850 m.resize(numberOfDofs, 1);
852 auto constexpr kNodesPerElement = TElement::kNodes;
853 auto const nQuadPts = eg.size();
854 for (
auto g = 0; g < nQuadPts; ++g)
856 auto const e = eg(g);
857 auto const nodes = E.col(e);
861 Meg.template block<kNodesPerElement, kNodesPerElement>(0, g * kNodesPerElement);
864 for (
auto i = 0; i < kNodesPerElement; ++i)
866 for (
auto j = 0; j < kNodesPerElement; ++j)
868 for (
auto d = 0; d < dims; ++d)
870 auto const ni = dims * nodes(i) + d;
FEM shape functions and gradients.
Eigen adaptors for ranges.
Finite Element Method (FEM)
Definition Concepts.h:19
auto ElementMassMatrices(Eigen::MatrixBase< TN > const &Neg, Eigen::DenseBase< TDerivedwg > const &wg, Eigen::DenseBase< TDerivedrhog > const &rhog)
Compute a matrix of horizontally stacked element mass matrices for all quadrature points.
Definition Mass.h:103
void ToElementMassMatrices(Eigen::MatrixBase< TN > const &Neg, Eigen::DenseBase< TDerivedwg > const &wg, Eigen::DenseBase< TDerivedrhog > const &rhog, Eigen::PlainObjectBase< TDerivedOut > &Meg)
Compute a matrix of horizontally stacked element mass matrices for all quadrature points.
Definition Mass.h:71
auto MassMatrix(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::DenseBase< TDerivedwg > const &wg, Eigen::DenseBase< TDerivedrhog > const &rhog, Eigen::MatrixBase< TDerivedNeg > const &Neg, int dims=1) -> Eigen::SparseMatrix< typename TDerivedNeg::Scalar, Options, typename TDerivedE::Scalar >
Construct the mass matrix operator's sparse matrix representation.
Definition Mass.h:664
auto LumpedMassMatrix(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::DenseBase< TDerivedwg > const &wg, Eigen::DenseBase< TDerivedrhog > const &rhog, Eigen::MatrixBase< TDerivedNeg > const &Neg, int dims=1) -> Eigen::Vector< typename TDerivedNeg::Scalar, Eigen::Dynamic >
Compute lumped mass matrix's diagonal vector (allocates output vector)
Definition Mass.h:500
auto ElementMassMatrix(Eigen::MatrixBase< TN > const &N, typename TN::Scalar w, typename TN::Scalar rho)
Compute the element mass matrix at a single quadrature point.
Definition Mass.h:43
void GemmMass(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::MatrixBase< TDerivedMe > const &Me, int dims, Eigen::MatrixBase< TDerivedIn > const &X, Eigen::DenseBase< TDerivedOut > &Y)
Compute mass matrix-matrix multiply using precomputed element mass matrices.
Definition Mass.h:597
void ToLumpedMassMatrix(Eigen::DenseBase< TDerivedE > const &E, Eigen::Index nNodes, Eigen::DenseBase< TDerivedeg > const &eg, Eigen::DenseBase< TDerivedwg > const &wg, Eigen::DenseBase< TDerivedrhog > const &rhog, Eigen::MatrixBase< TDerivedNeg > const &Neg, int dims, Eigen::PlainObjectBase< TDerivedOut > &m)
Compute lumped mass matrix's diagonal vector into existing output vector.
Definition Mass.h:791
Profiling utilities for the Physics-Based Animation Toolkit (PBAT)