Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62916157
MatrixPartitioned.hpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, May 16, 12:21
Size
10 KB
Mime Type
text/x-c++
Expires
Sat, May 18, 12:21 (2 d)
Engine
blob
Format
Raw Data
Handle
17707188
Attached To
rGOOSEFEM GooseFEM
MatrixPartitioned.hpp
View Options
/*
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
*/
#ifndef GOOSEFEM_MATRIXPARTITIONED_HPP
#define GOOSEFEM_MATRIXPARTITIONED_HPP
#include "MatrixPartitioned.h"
#include "Mesh.h"
namespace GooseFEM {
inline MatrixPartitioned::MatrixPartitioned(
const xt::xtensor<size_t, 2>& conn,
const xt::xtensor<size_t, 2>& dofs,
const xt::xtensor<size_t, 1>& iip)
: m_conn(conn), m_dofs(dofs), m_iip(iip)
{
m_nelem = m_conn.shape(0);
m_nne = m_conn.shape(1);
m_nnode = m_dofs.shape(0);
m_ndim = m_dofs.shape(1);
m_iiu = xt::setdiff1d(dofs, iip);
m_ndof = xt::amax(m_dofs)[0] + 1;
m_nnp = m_iip.size();
m_nnu = m_iiu.size();
m_part = Mesh::Reorder({m_iiu, m_iip}).get(m_dofs);
m_Tuu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
m_Tup.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
m_Tpu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
m_Tpp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
m_Auu.resize(m_nnu, m_nnu);
m_Aup.resize(m_nnu, m_nnp);
m_Apu.resize(m_nnp, m_nnu);
m_App.resize(m_nnp, m_nnp);
GOOSEFEM_ASSERT(xt::amax(m_conn)[0] + 1 == m_nnode);
GOOSEFEM_ASSERT(xt::amax(m_iip)[0] <= xt::amax(m_dofs)[0]);
GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
}
inline size_t MatrixPartitioned::nelem() const
{
return m_nelem;
}
inline size_t MatrixPartitioned::nne() const
{
return m_nne;
}
inline size_t MatrixPartitioned::nnode() const
{
return m_nnode;
}
inline size_t MatrixPartitioned::ndim() const
{
return m_ndim;
}
inline size_t MatrixPartitioned::ndof() const
{
return m_ndof;
}
inline size_t MatrixPartitioned::nnu() const
{
return m_nnu;
}
inline size_t MatrixPartitioned::nnp() const
{
return m_nnp;
}
inline xt::xtensor<size_t, 2> MatrixPartitioned::dofs() const
{
return m_dofs;
}
inline xt::xtensor<size_t, 1> MatrixPartitioned::iiu() const
{
return m_iiu;
}
inline xt::xtensor<size_t, 1> MatrixPartitioned::iip() const
{
return m_iip;
}
inline void MatrixPartitioned::assemble(const xt::xtensor<double, 3>& elemmat)
{
GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
m_Tuu.clear();
m_Tup.clear();
m_Tpu.clear();
m_Tpp.clear();
for (size_t e = 0; e < m_nelem; ++e) {
for (size_t m = 0; m < m_nne; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
size_t di = m_part(m_conn(e, m), i);
for (size_t n = 0; n < m_nne; ++n) {
for (size_t j = 0; j < m_ndim; ++j) {
size_t dj = m_part(m_conn(e, n), j);
if (di < m_nnu && dj < m_nnu) {
m_Tuu.push_back(Eigen::Triplet<double>(
di, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
}
else if (di < m_nnu) {
m_Tup.push_back(Eigen::Triplet<double>(
di, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
}
else if (dj < m_nnu) {
m_Tpu.push_back(Eigen::Triplet<double>(
di - m_nnu, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
}
else {
m_Tpp.push_back(Eigen::Triplet<double>(
di - m_nnu,
dj - m_nnu,
elemmat(e, m * m_ndim + i, n * m_ndim + j)));
}
}
}
}
}
}
m_Auu.setFromTriplets(m_Tuu.begin(), m_Tuu.end());
m_Aup.setFromTriplets(m_Tup.begin(), m_Tup.end());
m_Apu.setFromTriplets(m_Tpu.begin(), m_Tpu.end());
m_App.setFromTriplets(m_Tpp.begin(), m_Tpp.end());
m_changed = true;
}
inline void
MatrixPartitioned::reaction(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const
{
GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
Eigen::VectorXd X_u = this->asDofs_u(x);
Eigen::VectorXd X_p = this->asDofs_p(x);
Eigen::VectorXd B_p = m_Apu * X_u + m_App * X_p;
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) >= m_nnu) {
b(m, i) = B_p(m_part(m, i) - m_nnu);
}
}
}
}
inline void
MatrixPartitioned::reaction(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const
{
GOOSEFEM_ASSERT(x.size() == m_ndof);
GOOSEFEM_ASSERT(b.size() == m_ndof);
Eigen::VectorXd X_u = this->asDofs_u(x);
Eigen::VectorXd X_p = this->asDofs_p(x);
Eigen::VectorXd B_p = m_Apu * X_u + m_App * X_p;
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
b(m_iip(d)) = B_p(d);
}
}
inline void MatrixPartitioned::reaction_p(
const xt::xtensor<double, 1>& x_u,
const xt::xtensor<double, 1>& x_p,
xt::xtensor<double, 1>& b_p) const
{
GOOSEFEM_ASSERT(x_u.size() == m_nnu);
GOOSEFEM_ASSERT(x_p.size() == m_nnp);
GOOSEFEM_ASSERT(b_p.size() == m_nnp);
Eigen::VectorXd X_u(m_nnu, 1);
Eigen::VectorXd X_p(m_nnp, 1);
std::copy(x_u.begin(), x_u.end(), X_u.data());
std::copy(x_p.begin(), x_p.end(), X_p.data());
Eigen::VectorXd B_p = m_Apu * X_u + m_App * X_p;
std::copy(B_p.data(), B_p.data() + m_nnp, b_p.begin());
}
inline xt::xtensor<double, 2>
MatrixPartitioned::Reaction(const xt::xtensor<double, 2>& x, const xt::xtensor<double, 2>& b) const
{
xt::xtensor<double, 2> out = b;
this->reaction(x, out);
return out;
}
inline xt::xtensor<double, 1>
MatrixPartitioned::Reaction(const xt::xtensor<double, 1>& x, const xt::xtensor<double, 1>& b) const
{
xt::xtensor<double, 1> out = b;
this->reaction(x, out);
return out;
}
inline xt::xtensor<double, 1> MatrixPartitioned::Reaction_p(
const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const
{
xt::xtensor<double, 1> b_p = xt::empty<double>({m_nnp});
this->reaction_p(x_u, x_p, b_p);
return b_p;
}
inline Eigen::VectorXd MatrixPartitioned::asDofs_u(const xt::xtensor<double, 1>& dofval) const
{
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
Eigen::VectorXd dofval_u(m_nnu, 1);
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
dofval_u(d) = dofval(m_iiu(d));
}
return dofval_u;
}
inline Eigen::VectorXd MatrixPartitioned::asDofs_u(const xt::xtensor<double, 2>& nodevec) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
Eigen::VectorXd dofval_u(m_nnu, 1);
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) < m_nnu) {
dofval_u(m_part(m, i)) = nodevec(m, i);
}
}
}
return dofval_u;
}
inline Eigen::VectorXd MatrixPartitioned::asDofs_p(const xt::xtensor<double, 1>& dofval) const
{
GOOSEFEM_ASSERT(dofval.size() == m_ndof);
Eigen::VectorXd dofval_p(m_nnp, 1);
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
dofval_p(d) = dofval(m_iip(d));
}
return dofval_p;
}
inline Eigen::VectorXd MatrixPartitioned::asDofs_p(const xt::xtensor<double, 2>& nodevec) const
{
GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
Eigen::VectorXd dofval_p(m_nnp, 1);
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
if (m_part(m, i) >= m_nnu) {
dofval_p(m_part(m, i) - m_nnu) = nodevec(m, i);
}
}
}
return dofval_p;
}
template <class Solver>
inline void MatrixPartitionedSolver<Solver>::factorize(MatrixPartitioned& matrix)
{
if (!matrix.m_changed && !m_factor) {
return;
}
m_solver.compute(matrix.m_Auu);
m_factor = false;
matrix.m_changed = false;
}
template <class Solver>
inline void MatrixPartitionedSolver<Solver>::solve(
MatrixPartitioned& matrix, const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x)
{
GOOSEFEM_ASSERT(xt::has_shape(b, {matrix.m_nnode, matrix.m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(x, {matrix.m_nnode, matrix.m_ndim}));
this->factorize(matrix);
Eigen::VectorXd B_u = matrix.asDofs_u(b);
Eigen::VectorXd X_p = matrix.asDofs_p(x);
Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_Aup * X_p));
#pragma omp parallel for
for (size_t m = 0; m < matrix.m_nnode; ++m) {
for (size_t i = 0; i < matrix.m_ndim; ++i) {
if (matrix.m_part(m, i) < matrix.m_nnu) {
x(m, i) = X_u(matrix.m_part(m, i));
}
}
}
}
template <class Solver>
inline void MatrixPartitionedSolver<Solver>::solve(
MatrixPartitioned& matrix, const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x)
{
GOOSEFEM_ASSERT(b.size() == matrix.m_ndof);
GOOSEFEM_ASSERT(x.size() == matrix.m_ndof);
this->factorize(matrix);
Eigen::VectorXd B_u = matrix.asDofs_u(b);
Eigen::VectorXd X_p = matrix.asDofs_p(x);
Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_Aup * X_p));
#pragma omp parallel for
for (size_t d = 0; d < matrix.m_nnu; ++d) {
x(matrix.m_iiu(d)) = X_u(d);
}
}
template <class Solver>
inline void MatrixPartitionedSolver<Solver>::solve_u(
MatrixPartitioned& matrix,
const xt::xtensor<double, 1>& b_u,
const xt::xtensor<double, 1>& x_p,
xt::xtensor<double, 1>& x_u)
{
GOOSEFEM_ASSERT(b_u.size() == matrix.m_nnu);
GOOSEFEM_ASSERT(x_p.size() == matrix.m_nnp);
GOOSEFEM_ASSERT(x_u.size() == matrix.m_nnu);
this->factorize(matrix);
Eigen::VectorXd B_u(matrix.m_nnu, 1);
Eigen::VectorXd X_p(matrix.m_nnp, 1);
std::copy(b_u.begin(), b_u.end(), B_u.data());
std::copy(x_p.begin(), x_p.end(), X_p.data());
Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_Aup * X_p));
std::copy(X_u.data(), X_u.data() + matrix.m_nnu, x_u.begin());
}
template <class Solver>
inline xt::xtensor<double, 2> MatrixPartitionedSolver<Solver>::Solve(
MatrixPartitioned& matrix, const xt::xtensor<double, 2>& b, const xt::xtensor<double, 2>& x)
{
xt::xtensor<double, 2> out = x;
this->solve(matrix, b, out);
return out;
}
template <class Solver>
inline xt::xtensor<double, 1> MatrixPartitionedSolver<Solver>::Solve(
MatrixPartitioned& matrix, const xt::xtensor<double, 1>& b, const xt::xtensor<double, 1>& x)
{
xt::xtensor<double, 1> out = x;
this->solve(matrix, b, out);
return out;
}
template <class Solver>
inline xt::xtensor<double, 1> MatrixPartitionedSolver<Solver>::Solve_u(
MatrixPartitioned& matrix, const xt::xtensor<double, 1>& b_u, const xt::xtensor<double, 1>& x_p)
{
xt::xtensor<double, 1> x_u = xt::empty<double>({matrix.m_nnu});
this->solve_u(matrix, b_u, x_p, x_u);
return x_u;
}
} // namespace GooseFEM
#endif
Event Timeline
Log In to Comment