Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F59579551
MatrixDiagonalPartitioned.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
Tue, Apr 23, 16:07
Size
8 KB
Mime Type
text/x-c
Expires
Thu, Apr 25, 16:07 (2 d)
Engine
blob
Format
Raw Data
Handle
17213409
Attached To
rGOOSEFEM GooseFEM
MatrixDiagonalPartitioned.hpp
View Options
/**
Implementation of MatrixDiagonalPartitioned.h
\file MatrixDiagonalPartitioned.hpp
\copyright Copyright 2017. Tom de Geus. All rights reserved.
\license This project is released under the GNU Public License (GPLv3).
*/
#ifndef GOOSEFEM_MATRIXDIAGONALPARTITIONED_HPP
#define GOOSEFEM_MATRIXDIAGONALPARTITIONED_HPP
#include "MatrixDiagonalPartitioned.h"
#include "Mesh.h"
namespace GooseFEM {
inline MatrixDiagonalPartitioned::MatrixDiagonalPartitioned(
const xt::xtensor<size_t, 2>& conn,
const xt::xtensor<size_t, 2>& dofs,
const xt::xtensor<size_t, 1>& iip)
: MatrixDiagonal(conn, dofs)
{
m_iip = iip;
m_iiu = xt::setdiff1d(dofs, iip);
m_nnp = m_iip.size();
m_nnu = m_iiu.size();
m_part = Mesh::Reorder({m_iiu, m_iip}).apply(m_dofs);
m_Auu = xt::empty<double>({m_nnu});
m_App = xt::empty<double>({m_nnp});
m_inv_uu = xt::empty<double>({m_nnu});
GOOSEFEM_ASSERT(xt::amax(m_iip)() <= xt::amax(m_dofs)());
}
inline size_t MatrixDiagonalPartitioned::nnu() const
{
return m_nnu;
}
inline size_t MatrixDiagonalPartitioned::nnp() const
{
return m_nnp;
}
inline xt::xtensor<size_t, 1> MatrixDiagonalPartitioned::iiu() const
{
return m_iiu;
}
inline xt::xtensor<size_t, 1> MatrixDiagonalPartitioned::iip() const
{
return m_iip;
}
inline void MatrixDiagonalPartitioned::factorize()
{
if (!m_changed) {
return;
}
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
m_inv_uu(d) = 1.0 / m_Auu(d);
}
m_changed = false;
}
inline void MatrixDiagonalPartitioned::assemble(const xt::xtensor<double, 3>& elemmat)
{
GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
GOOSEFEM_ASSERT(Element::isDiagonal(elemmat));
m_Auu.fill(0.0);
m_App.fill(0.0);
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 d = m_part(m_conn(e, m), i);
if (d < m_nnu) {
m_Auu(d) += elemmat(e, m * m_ndim + i, m * m_ndim + i);
}
else {
m_App(d - m_nnu) += elemmat(e, m * m_ndim + i, m * m_ndim + i);
}
}
}
}
m_changed = true;
}
inline void MatrixDiagonalPartitioned::set(const xt::xtensor<double, 1>& A)
{
GOOSEFEM_ASSERT(A.size() == m_ndof);
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
m_Auu(d) = A(m_iiu(d));
}
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
m_App(d) = A(m_iip(d));
}
m_changed = true;
}
inline void
MatrixDiagonalPartitioned::dot(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}));
#pragma omp parallel for
for (size_t m = 0; m < m_nnode; ++m) {
for (size_t i = 0; i < m_ndim; ++i) {
size_t d = m_part(m, i);
if (d < m_nnu) {
b(m, i) = m_Auu(d) * x(m, i);
}
else {
b(m, i) = m_App(d - m_nnu) * x(m, i);
}
}
}
}
inline void
MatrixDiagonalPartitioned::dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const
{
GOOSEFEM_ASSERT(x.size() == m_ndof);
GOOSEFEM_ASSERT(b.size() == m_ndof);
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
b(m_iiu(d)) = m_Auu(d) * x(m_iiu(d));
}
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
b(m_iip(d)) = m_App(d) * x(m_iip(d));
}
}
inline void MatrixDiagonalPartitioned::dot_u(
const xt::xtensor<double, 1>& x_u,
const xt::xtensor<double, 1>& x_p,
xt::xtensor<double, 1>& b_u) const
{
UNUSED(x_p);
GOOSEFEM_ASSERT(x_u.size() == m_nnu);
GOOSEFEM_ASSERT(x_p.size() == m_nnp);
GOOSEFEM_ASSERT(b_u.size() == m_nnu);
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
b_u(d) = m_Auu(d) * x_u(d);
}
}
inline void MatrixDiagonalPartitioned::dot_p(
const xt::xtensor<double, 1>& x_u,
const xt::xtensor<double, 1>& x_p,
xt::xtensor<double, 1>& b_p) const
{
UNUSED(x_u);
GOOSEFEM_ASSERT(x_u.size() == m_nnu);
GOOSEFEM_ASSERT(x_p.size() == m_nnp);
GOOSEFEM_ASSERT(b_p.size() == m_nnp);
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
b_p(d) = m_App(d) * x_p(d);
}
}
inline void
MatrixDiagonalPartitioned::solve(const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x)
{
GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
this->factorize();
#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) {
x(m, i) = m_inv_uu(m_part(m, i)) * b(m, i);
}
}
}
}
inline void
MatrixDiagonalPartitioned::solve(const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x)
{
GOOSEFEM_ASSERT(b.size() == m_ndof);
GOOSEFEM_ASSERT(x.size() == m_ndof);
this->factorize();
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
x(m_iiu(d)) = m_inv_uu(d) * b(m_iiu(d));
}
}
inline void MatrixDiagonalPartitioned::solve_u(
const xt::xtensor<double, 1>& b_u,
const xt::xtensor<double, 1>& x_p,
xt::xtensor<double, 1>& x_u)
{
UNUSED(x_p);
GOOSEFEM_ASSERT(b_u.size() == m_nnu);
GOOSEFEM_ASSERT(x_p.size() == m_nnp);
GOOSEFEM_ASSERT(x_u.size() == m_nnu);
this->factorize();
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
x_u(d) = m_inv_uu(d) * b_u(d);
}
}
inline void MatrixDiagonalPartitioned::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}));
#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) = m_App(m_part(m, i) - m_nnu) * x(m, i);
}
}
}
}
inline void MatrixDiagonalPartitioned::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);
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
b(m_iip(d)) = m_App(d) * x(m_iip(d));
}
}
inline void MatrixDiagonalPartitioned::reaction_p(
const xt::xtensor<double, 1>& x_u,
const xt::xtensor<double, 1>& x_p,
xt::xtensor<double, 1>& b_p) const
{
UNUSED(x_u);
GOOSEFEM_ASSERT(x_u.size() == m_nnu);
GOOSEFEM_ASSERT(x_p.size() == m_nnp);
GOOSEFEM_ASSERT(b_p.size() == m_nnp);
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
b_p(d) = m_App(d) * x_p(d);
}
}
inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Todiagonal() const
{
xt::xtensor<double, 1> ret = xt::zeros<double>({m_ndof});
#pragma omp parallel for
for (size_t d = 0; d < m_nnu; ++d) {
ret(m_iiu(d)) = m_Auu(d);
}
#pragma omp parallel for
for (size_t d = 0; d < m_nnp; ++d) {
ret(m_iip(d)) = m_App(d);
}
return ret;
}
inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Dot_u(
const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const
{
xt::xtensor<double, 1> b_u = xt::empty<double>({m_nnu});
this->dot_u(x_u, x_p, b_u);
return b_u;
}
inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Dot_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->dot_p(x_u, x_p, b_p);
return b_p;
}
inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Solve_u(
const xt::xtensor<double, 1>& b_u, const xt::xtensor<double, 1>& x_p)
{
xt::xtensor<double, 1> x_u = xt::empty<double>({m_nnu});
this->solve_u(b_u, x_p, x_u);
return x_u;
}
inline xt::xtensor<double, 2> MatrixDiagonalPartitioned::Reaction(
const xt::xtensor<double, 2>& x, const xt::xtensor<double, 2>& b) const
{
xt::xtensor<double, 2> ret = b;
this->reaction(x, ret);
return ret;
}
inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Reaction(
const xt::xtensor<double, 1>& x, const xt::xtensor<double, 1>& b) const
{
xt::xtensor<double, 1> ret = b;
this->reaction(x, ret);
return ret;
}
inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::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;
}
} // namespace GooseFEM
#endif
Event Timeline
Log In to Comment