1713 {
1715 if (row_type != MBTET)
1717
1718 if (auto ep_fe_vol_ptr =
1719 dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
1720 getFEMethod())) {
1721
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
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
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()) {
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
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()) {
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) {
1800 &*cind.begin(), &*
m.data().begin(), ADD_VALUES);
1801
1803 };
1804
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);
1816 const auto info =
lapack_dsysv(
'L', nb, nb, &*
m.data().begin(), nb,
1817 &*
iPIV.begin(), &*inv.data().begin(), nb,
1819
1820
1821
1822
1824 };
1825
1827 std::string field, auto &inv) {
1829
1831 bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
1832 if (
bit != bmc.get<1>().end()) {
1833
1835 CHKERR invert_symm_mat(
m, inv);
1836
1837 } else
1839 "%s matrix not found", field.c_str());
1840
1842 };
1843
1846
1847 const int nb =
m.size1();
1848
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 =
1859 &*trans_inv.data().begin(), nb);
1860 if (info != 0)
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 =
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
1880 CHKERR invert_nonsymm_mat(
m, inv);
1881
1883 std::cerr << prod(
m, inv) << endl;
1884 std::cerr << endl;
1885 }
1886
1887 } else
1889 "%s matrix not found", field.c_str());
1890
1892 };
1893
1894 auto create_block_schur =
1899
1900 for (
auto &
bit : add_bmc) {
1901 bit.unSetAtElement();
1903 }
1904
1905 for (
auto &
bit : bmc) {
1906 if (
bit.setAtElement &&
bit.rowField != field &&
1907 bit.colField != field) {
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,
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) {
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) {
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 =
1954 bool debug =
false) {
1956 for (
auto &
bit : add_bmc) {
1957 if (
bit.setAtElement)
1959 }
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 =
1979
1980 for (
auto &
bit : add_bmc) {
1981 bit.unSetAtElement();
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,
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,
2001 bit->rowInd,
bit->colInd);
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);
2012 break;
2013 }
2014 }
2015 if (row_it == bmc.get<3>().end())
2017 "row field not found %s", field.c_str());
2018 }
2019
2021 };
2022
2028
2029 if (SBB) {
2036
2037 if (SOO) {
2045 } else {
2047 }
2049
2050 if (Sww) {
2056 }
2057 }
2058 }
2059 }
2060 }
2061
2063}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'm', SPACE_DIM > m
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)
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)
UBlasVector< int > VectorInt
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
MatrixDouble K
local tangent matrix
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