1715 if (row_type != MBTET)
1718 if (
auto ep_fe_vol_ptr =
1719 dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *
>(
1732 auto set_data = [&](
auto fe) {
1734 aoSBB = fe->aoSBubble;
1735 aoSOO = fe->aoSOmega;
1743 set_data(ep_fe_vol_ptr);
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));
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();
1767 auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
1768 row_name, col_name, row_type, col_type, row_side, col_side,
m,
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();
1787 auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
1788 row_name, col_name, row_type, col_type, row_side, col_side,
m,
1794 auto assemble_block = [&](
auto &bit, Mat S) {
1800 &*cind.begin(), &*
m.data().begin(), ADD_VALUES);
1807 const int nb =
m.size1();
1809 inv.resize(nb, nb,
false);
1811 for (
int cc = 0; cc != nb; ++cc)
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,
1827 std::string field,
auto &inv) {
1831 bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
1832 if (bit != bmc.get<1>().end()) {
1835 CHKERR invert_symm_mat(
m, inv);
1839 "%s matrix not found", field.c_str());
1847 const int nb =
m.size1();
1851 trans_inv.resize(nb, nb,
false);
1853 for (
int c = 0;
c != nb; ++
c)
1854 trans_inv(
c,
c) = -1;
1856 iPIV.resize(nb,
false);
1859 &*trans_inv.data().begin(), nb);
1862 "Can not invert matrix info = %d", info);
1864 inv.resize(nb, nb,
false);
1865 noalias(inv) = trans(trans_inv);
1870 auto invert_nonsymm_schur =
1872 auto &inv,
const bool debug =
false) {
1875 auto bit = bmc.get<1>().find(
1876 boost::make_tuple(field, field, MBTET, MBTET));
1877 if (bit != bmc.get<1>().end()) {
1880 CHKERR invert_nonsymm_mat(
m, inv);
1883 std::cerr << prod(
m, inv) << endl;
1889 "%s matrix not found", field.c_str());
1894 auto create_block_schur =
1900 for (
auto &bit : add_bmc) {
1901 bit.unSetAtElement();
1905 for (
auto &bit : bmc) {
1906 if (bit.setAtElement && bit.rowField != field &&
1907 bit.colField != field) {
1912 CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
1913 CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1915 auto it = set_block(add_bmc, bit.rowField, bit.colField,
1916 bit.rowSide, bit.colSide, bit.rowType,
1917 bit.colType,
m, rind, cind);
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);
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);
1938 CHKERR AOApplicationToPetsc(ao, rind.size(),
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);
1952 auto assemble_schur =
1954 bool debug =
false) {
1956 for (
auto &bit : add_bmc) {
1957 if (bit.setAtElement)
1958 CHKERR assemble_block(bit, S);
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;
1968 std::cerr << std::endl;
1973 auto precondition_schur =
1980 for (
auto &bit : add_bmc) {
1981 bit.unSetAtElement();
1985 for (
auto &bit : bmc) {
1986 if (bit.setAtElement) {
1987 if (bit.rowField != field || bit.colField != field)
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);
1995 auto bit = bmc.get<1>().find(
1996 boost::make_tuple(field, field, MBTET, MBTET));
1997 if (bit->setAtElement && bit != bmc.get<1>().end()) {
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);
2003 m +=
eps * diag_mat;
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);
2015 if (row_it == bmc.get<3>().end())
2017 "row field not found %s", field.c_str());
Public Member Functions inherited from EshelbianPlasticity::OpAssembleVolume