v0.13.2
Loading...
Searching...
No Matches
Public Member Functions | Private Attributes | List of all members
EshelbianPlasticity::OpSpatialSchurEnd Struct Reference

#include <users_modules/eshelbian_plasticty/src/EshelbianPlasticity.hpp>

Inheritance diagram for EshelbianPlasticity::OpSpatialSchurEnd:
[legend]
Collaboration diagram for EshelbianPlasticity::OpSpatialSchurEnd:
[legend]

Public Member Functions

 OpSpatialSchurEnd (const std::string &row_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const double eps)
 
MoFEMErrorCode integrate (EntData &row_data)
 
MoFEMErrorCode assemble (int side, EntityType type, EntData &data)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 
- Public Member Functions inherited from EshelbianPlasticity::OpAssembleVolume
 OpAssembleVolume (const std::string &field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type)
 
 OpAssembleVolume (const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type, const bool assemble_symmetry)
 
- Public Member Functions inherited from EshelbianPlasticity::OpAssembleBasic< VolUserDataOperator >
 OpAssembleBasic (const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type)
 
 OpAssembleBasic (const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type, const bool assemble_symmetry)
 
virtual MoFEMErrorCode integrate (EntData &data)
 
virtual MoFEMErrorCode integrate (int row_side, EntityType row_type, EntData &data)
 
virtual MoFEMErrorCode integrate (EntData &row_data, EntData &col_data)
 
virtual MoFEMErrorCode assemble (EntData &data)
 
virtual MoFEMErrorCode assemble (int row_side, EntityType row_type, EntData &data)
 
virtual MoFEMErrorCode assemble (int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 
MoFEMErrorCode doWork (int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
 

Private Attributes

const double eps
 
MatrixDouble invMat
 
VectorInt iPIV
 
VectorDouble lapackWork
 
map< std::string, MatrixDouble > invBlockMat
 
map< std::string, DataAtIntegrationPts::BlockMatContainorblockMat
 

Additional Inherited Members

- Public Attributes inherited from EshelbianPlasticity::OpAssembleBasic< VolUserDataOperator >
const bool assembleSymmetry
 
boost::shared_ptr< DataAtIntegrationPtsdataAtPts
 data at integration pts More...
 
VectorDouble nF
 local right hand side vector More...
 
MatrixDouble K
 local tangent matrix More...
 
MatrixDouble transposeK
 

Detailed Description

Definition at line 1108 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

◆ OpSpatialSchurEnd()

EshelbianPlasticity::OpSpatialSchurEnd::OpSpatialSchurEnd ( const std::string &  row_field,
boost::shared_ptr< DataAtIntegrationPts > &  data_ptr,
const double  eps 
)
inline

Definition at line 1109 of file EshelbianPlasticity.hpp.

1112 : OpAssembleVolume(row_field, data_ptr, OPROW), eps(eps) {
1113 sYmm = false;
1114 }
OpAssembleVolume(const std::string &field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type)

Member Function Documentation

◆ assemble()

MoFEMErrorCode EshelbianPlasticity::OpSpatialSchurEnd::assemble ( int  side,
EntityType  type,
EntData data 
)
virtual

Reimplemented from EshelbianPlasticity::OpAssembleBasic< VolUserDataOperator >.

Examples
EshelbianOperators.cpp.

Definition at line 1712 of file EshelbianOperators.cpp.

1713 {
1715 if (row_type != MBTET)
1717
1718 if (auto ep_fe_vol_ptr =
1719 dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
1720 getFEMethod())) {
1721
1722 SmartPetscObj<AO> aoSuu;
1723 SmartPetscObj<AO> aoSBB;
1724 SmartPetscObj<AO> aoSOO;
1725 SmartPetscObj<AO> aoSww;
1726
1731
1732 auto set_data = [&](auto fe) {
1733 aoSuu = fe->aoSuu;
1734 aoSBB = fe->aoSBubble;
1735 aoSOO = fe->aoSOmega;
1736 aoSww = fe->aoSw;
1737
1738 Suu = fe->Suu;
1739 SBB = fe->SBubble;
1740 SOO = fe->SOmega;
1741 Sww = fe->Sw;
1742 };
1743 set_data(ep_fe_vol_ptr);
1744
1745 if (Suu) {
1746
1747 auto find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1748 auto &row_name, auto &col_name, auto row_side,
1749 auto col_side, auto row_type, auto col_type) {
1750 return add_bmc.get<0>().find(boost::make_tuple(
1751 row_name, col_name, row_type, col_type, row_side, col_side));
1752 };
1753
1754 auto set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1755 auto &row_name, auto &col_name, auto row_side,
1756 auto col_side, auto row_type, auto col_type,
1757 const auto &m, const auto &row_ind,
1758 const auto &col_ind) {
1759 auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
1760 row_type, col_type);
1761 if (it != add_bmc.get<0>().end()) {
1762 it->setMat(m);
1763 it->setInd(row_ind, col_ind);
1764 it->setSetAtElement();
1765 return it;
1766 } else {
1767 auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
1768 row_name, col_name, row_type, col_type, row_side, col_side, m,
1769 row_ind, col_ind));
1770 return p_it.first;
1771 }
1772 };
1773
1774 auto add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1775 auto &row_name, auto &col_name, auto row_side,
1776 auto col_side, auto row_type, auto col_type,
1777 const auto &m, const auto &row_ind,
1778 const auto &col_ind) {
1779 auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
1780 row_type, col_type);
1781 if (it != add_bmc.get<0>().end()) {
1782 it->addMat(m);
1783 it->setInd(row_ind, col_ind);
1784 it->setSetAtElement();
1785 return it;
1786 } else {
1787 auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
1788 row_name, col_name, row_type, col_type, row_side, col_side, m,
1789 row_ind, col_ind));
1790 return p_it.first;
1791 }
1792 };
1793
1794 auto assemble_block = [&](auto &bit, Mat S) {
1796 const VectorInt &rind = bit.rowInd;
1797 const VectorInt &cind = bit.colInd;
1798 const MatrixDouble &m = bit.M;
1799 CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
1800 &*cind.begin(), &*m.data().begin(), ADD_VALUES);
1801
1803 };
1804
1805 auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
1807 const int nb = m.size1();
1808
1809 inv.resize(nb, nb, false);
1810 inv.clear();
1811 for (int cc = 0; cc != nb; ++cc)
1812 inv(cc, cc) = -1;
1813
1814 iPIV.resize(nb, false);
1815 lapackWork.resize(nb * nb, false);
1816 const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
1817 &*iPIV.begin(), &*inv.data().begin(), nb,
1818 &*lapackWork.begin(), nb * nb);
1819 // if (info != 0)
1820 // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1821 // "Can not invert matrix info = %d", info);
1822
1824 };
1825
1826 auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
1827 std::string field, auto &inv) {
1829
1830 auto bit =
1831 bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
1832 if (bit != bmc.get<1>().end()) {
1833
1834 auto &m = *const_cast<MatrixDouble *>(&(bit->M));
1835 CHKERR invert_symm_mat(m, inv);
1836
1837 } else
1838 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1839 "%s matrix not found", field.c_str());
1840
1842 };
1843
1844 auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
1846
1847 const int nb = m.size1();
1848
1849 MatrixDouble trans_m = trans(m);
1850 MatrixDouble trans_inv;
1851 trans_inv.resize(nb, nb, false);
1852 trans_inv.clear();
1853 for (int c = 0; c != nb; ++c)
1854 trans_inv(c, c) = -1;
1855
1856 iPIV.resize(nb, false);
1857 const auto info =
1858 lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
1859 &*trans_inv.data().begin(), nb);
1860 if (info != 0)
1861 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1862 "Can not invert matrix info = %d", info);
1863
1864 inv.resize(nb, nb, false);
1865 noalias(inv) = trans(trans_inv);
1866
1868 };
1869
1870 auto invert_nonsymm_schur =
1871 [&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
1872 auto &inv, const bool debug = false) {
1874
1875 auto bit = bmc.get<1>().find(
1876 boost::make_tuple(field, field, MBTET, MBTET));
1877 if (bit != bmc.get<1>().end()) {
1878
1879 auto &m = *const_cast<MatrixDouble *>(&(bit->M));
1880 CHKERR invert_nonsymm_mat(m, inv);
1881
1882 if (debug) {
1883 std::cerr << prod(m, inv) << endl;
1884 std::cerr << endl;
1885 }
1886
1887 } else
1888 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1889 "%s matrix not found", field.c_str());
1890
1892 };
1893
1894 auto create_block_schur =
1897 std::string field, AO ao, MatrixDouble &inv) {
1899
1900 for (auto &bit : add_bmc) {
1901 bit.unSetAtElement();
1902 bit.clearMat();
1903 }
1904
1905 for (auto &bit : bmc) {
1906 if (bit.setAtElement && bit.rowField != field &&
1907 bit.colField != field) {
1908 VectorInt rind = bit.rowInd;
1909 VectorInt cind = bit.colInd;
1910 const MatrixDouble &m = bit.M;
1911 if (ao) {
1912 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
1913 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1914 }
1915 auto it = set_block(add_bmc, bit.rowField, bit.colField,
1916 bit.rowSide, bit.colSide, bit.rowType,
1917 bit.colType, m, rind, cind);
1918 }
1919 }
1920
1921 for (auto &bit_col : bmc) {
1922 if (bit_col.setAtElement && bit_col.rowField == field &&
1923 bit_col.colField != field) {
1924 const MatrixDouble &cm = bit_col.M;
1925 VectorInt cind = bit_col.colInd;
1926 invMat.resize(inv.size1(), cm.size2(), false);
1927 noalias(invMat) = prod(inv, cm);
1928 if (ao)
1929 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1930 for (auto &bit_row : bmc) {
1931 if (bit_row.setAtElement && bit_row.rowField != field &&
1932 bit_row.colField == field) {
1933 const MatrixDouble &rm = bit_row.M;
1934 VectorInt rind = bit_row.rowInd;
1935 K.resize(rind.size(), cind.size(), false);
1936 noalias(K) = prod(rm, invMat);
1937 if (ao)
1938 CHKERR AOApplicationToPetsc(ao, rind.size(),
1939 &*rind.begin());
1940 auto it = add_block(add_bmc, bit_row.rowField,
1941 bit_col.colField, bit_row.rowSide,
1942 bit_col.colSide, bit_row.rowType,
1943 bit_col.colType, K, rind, cind);
1944 }
1945 }
1946 }
1947 }
1948
1950 };
1951
1952 auto assemble_schur =
1953 [&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
1954 bool debug = false) {
1956 for (auto &bit : add_bmc) {
1957 if (bit.setAtElement)
1958 CHKERR assemble_block(bit, S);
1959 }
1960 if (debug) {
1961 for (auto &bit : add_bmc) {
1962 if (bit.setAtElement) {
1963 std::cerr << "assemble: " << bit.rowField << " "
1964 << bit.colField << endl;
1965 std::cerr << bit.M << endl;
1966 }
1967 }
1968 std::cerr << std::endl;
1969 }
1971 };
1972
1973 auto precondition_schur =
1976 const std::string field, const MatrixDouble &diag_mat,
1977 const double eps) {
1979
1980 for (auto &bit : add_bmc) {
1981 bit.unSetAtElement();
1982 bit.clearMat();
1983 }
1984
1985 for (auto &bit : bmc) {
1986 if (bit.setAtElement) {
1987 if (bit.rowField != field || bit.colField != field)
1988 auto it =
1989 set_block(add_bmc, bit.rowField, bit.colField,
1990 bit.rowSide, bit.colSide, bit.rowType,
1991 bit.colType, bit.M, bit.rowInd, bit.colInd);
1992 }
1993 }
1994
1995 auto bit = bmc.get<1>().find(
1996 boost::make_tuple(field, field, MBTET, MBTET));
1997 if (bit->setAtElement && bit != bmc.get<1>().end()) {
1998 auto it =
1999 set_block(add_bmc, bit->rowField, bit->colField, bit->rowSide,
2000 bit->colSide, bit->rowType, bit->colType, bit->M,
2001 bit->rowInd, bit->colInd);
2002 MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2003 m += eps * diag_mat;
2004 } else {
2005 auto row_it = bmc.get<3>().lower_bound(field);
2006 for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
2007 if (row_it->setAtElement) {
2008 auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
2009 diag_mat, row_it->rowInd, row_it->rowInd);
2010 MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2011 m *= eps;
2012 break;
2013 }
2014 }
2015 if (row_it == bmc.get<3>().end())
2016 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2017 "row field not found %s", field.c_str());
2018 }
2019
2021 };
2022
2023 CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
2024 invBlockMat["uu"]);
2025 CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
2026 "u", aoSuu, invBlockMat["uu"]);
2027 CHKERR assemble_schur(blockMat["uu"], Suu);
2028
2029 if (SBB) {
2030 CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
2031 CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
2032 aoSBB, invBlockMat["BB"]);
2033 CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
2034 *dataAtPts->ooMatPtr, eps);
2035 CHKERR assemble_schur(blockMat["precBBOO"], SBB);
2036
2037 if (SOO) {
2038 CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
2039 invBlockMat["OO"]);
2040 CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
2041 "omega", aoSOO, invBlockMat["OO"]);
2042 if (dataAtPts->wwMatPtr) {
2043 CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
2044 *dataAtPts->wwMatPtr, -eps);
2045 } else {
2046 blockMat["precOOww"] = blockMat["OO"];
2047 }
2048 CHKERR assemble_schur(blockMat["precOOww"], SOO);
2049
2050 if (Sww) {
2051 CHKERR invert_symm_schur(blockMat["precOOww"], "w",
2052 invBlockMat["ww"]);
2053 CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
2054 aoSww, invBlockMat["ww"]);
2055 CHKERR assemble_schur(blockMat["ww"], Sww);
2056 }
2057 }
2058 }
2059 }
2060 }
2061
2063}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
static const bool debug
FTensor::Index< 'm', SPACE_DIM > m
auto bit
set bit
const double c
speed of light (cm/ns)
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:176
static __CLPK_integer lapack_dsysv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:220
UBlasVector< int > VectorInt
Definition: Types.hpp:67
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
multi_index_container< BlockMatData, indexed_by< ordered_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType >, member< BlockMatData, int, &BlockMatData::rowSide >, member< BlockMatData, int, &BlockMatData::colSide > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField > > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::rowField > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::colField > > > > BlockMatContainor
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
map< std::string, MatrixDouble > invBlockMat
map< std::string, DataAtIntegrationPts::BlockMatContainor > blockMat
intrusive_ptr for managing petsc objects

◆ doWork()

MoFEMErrorCode EshelbianPlasticity::OpSpatialSchurEnd::doWork ( int  side,
EntityType  type,
EntData data 
)
inline

Definition at line 1118 of file EshelbianPlasticity.hpp.

1118 {
1119 return assemble(side, type, data);
1120 }
MoFEMErrorCode assemble(int side, EntityType type, EntData &data)

◆ integrate()

MoFEMErrorCode EshelbianPlasticity::OpSpatialSchurEnd::integrate ( EntData row_data)
inlinevirtual

Reimplemented from EshelbianPlasticity::OpAssembleBasic< VolUserDataOperator >.

Definition at line 1116 of file EshelbianPlasticity.hpp.

1116{ return 0; }

Member Data Documentation

◆ blockMat

map<std::string, DataAtIntegrationPts::BlockMatContainor> EshelbianPlasticity::OpSpatialSchurEnd::blockMat
private
Examples
EshelbianOperators.cpp.

Definition at line 1130 of file EshelbianPlasticity.hpp.

◆ eps

const double EshelbianPlasticity::OpSpatialSchurEnd::eps
private
Examples
EshelbianOperators.cpp.

Definition at line 1123 of file EshelbianPlasticity.hpp.

◆ invBlockMat

map<std::string, MatrixDouble> EshelbianPlasticity::OpSpatialSchurEnd::invBlockMat
private
Examples
EshelbianOperators.cpp.

Definition at line 1129 of file EshelbianPlasticity.hpp.

◆ invMat

MatrixDouble EshelbianPlasticity::OpSpatialSchurEnd::invMat
private
Examples
EshelbianOperators.cpp.

Definition at line 1125 of file EshelbianPlasticity.hpp.

◆ iPIV

VectorInt EshelbianPlasticity::OpSpatialSchurEnd::iPIV
private
Examples
EshelbianOperators.cpp.

Definition at line 1126 of file EshelbianPlasticity.hpp.

◆ lapackWork

VectorDouble EshelbianPlasticity::OpSpatialSchurEnd::lapackWork
private
Examples
EshelbianOperators.cpp.

Definition at line 1127 of file EshelbianPlasticity.hpp.


The documentation for this struct was generated from the following files: